專利名稱：用于在模具填充過程的模擬中描述粒子統(tǒng)計取向分布的方法和裝置的制作方法
技術(shù)領(lǐng)域：
本發(fā)明涉及對進入模具型腔的包含粒子的懸浮液的流進行三維建模的領(lǐng)域，并且更具體地，涉及用于在過程的模擬中描述非球形粒子的統(tǒng)計取向分布的方法和裝置，在該過程中用包含大量非球形粒子的懸浮液填充模具型腔。
背景技術(shù)：
注塑成型過程或者金屬鑄造過程的真實3-D模擬涉及多達幾十萬個方程的系統(tǒng)。在過去已經(jīng)取得了進展來改進模擬方法的效率以處理這些復(fù)雜的計算。使用優(yōu)化的軟件以及現(xiàn)代工作站的處理能力，可以在工作場所執(zhí)行這種模擬，即以足夠快的速度獲取結(jié)果以適應(yīng)純科學(xué)研究領(lǐng)域之外的領(lǐng)域，并且可以由注塑成型物件的研究和開發(fā)部門、鑄造廠和制造廠中的工程師應(yīng)用。
當(dāng)將粒子(諸如纖維)被加入到聚合物合成物，并且需要描述纖維的取向分布時，模擬和相關(guān)方程組明顯變得更加復(fù)雜。由于計算時間太長或者模擬的精度不夠，由于模擬的復(fù)雜性不允許在現(xiàn)今的工作站得到可接受的結(jié)果，因此迄今還沒有在工作場所成功地引入這種過程的真實3-D模擬。
在纖維增強部件中，通常對于開發(fā)工程師來說，具有纖維取向分布的描述以便能夠預(yù)測該部件的張力和翹曲方面是至關(guān)重要的。通常纖維用于改進塑料部件的機械性能。但是因此(熱)機械性能(如熱膨脹、強度和剛度)取決于纖維的取向。
近年來在許多行業(yè)中注塑成型塑料部件的使用已經(jīng)穩(wěn)定地增長。與以前相比，電子設(shè)備、消費品、醫(yī)療設(shè)備和汽車零件的制造商正用塑料制造越來越多他們的產(chǎn)品和在他們的產(chǎn)品中使用的部件。
因為注塑成型纖維增強部件提供提高的強度/重量比、耐用性、部件集成度和較低的成本，因此它們正在取代結(jié)構(gòu)金屬部件。
同時，競爭壓力正驅(qū)使塑料注塑成型工業(yè)中的制造商去尋找優(yōu)化設(shè)計的新方法以更好地使設(shè)計匹配制造過程。當(dāng)在設(shè)計開發(fā)過程中較晚發(fā)現(xiàn)了需要修改部件或者模具構(gòu)造時，用于實施必要的改變的延遲和相關(guān)成本明顯高于在設(shè)計開發(fā)階段的較早階段。
想保證他們的部件是可生產(chǎn)的并且將最優(yōu)地執(zhí)行的公司想要使用計算機輔助工程技術(shù)來對注塑模具中的復(fù)雜流和所得到的纖維取向進行模擬和建模，以便在設(shè)計階段早期更好地理解制造過程并且將這種認(rèn)知結(jié)合到部件設(shè)計中。
當(dāng)設(shè)計用于纖維增強部件的注塑模具以及要在其中制造的纖維增強部件時，應(yīng)該考慮若干因素。諸如整體部件幾何形狀、最小和最大壁厚度、模具中的通過其注入液態(tài)聚合物和纖維懸浮液的門的數(shù)量和位置、模具中的型腔中的氣體通過其排出的排氣口的數(shù)量和位置、聚合物成分和性質(zhì)、纖維性質(zhì)和數(shù)量、收縮、容差和纖維取向分布的參數(shù)是其中的一些。由于緊密相互關(guān)聯(lián)的關(guān)系，部件和模具設(shè)計不能可靠地純粹基于終端部件的形式和功能，而是應(yīng)該也考慮制造過程的影響。
計算機輔助工程模擬可以有利地用來為設(shè)計和制造工程師提供關(guān)于在注塑成型過程期間在模具型腔內(nèi)可能發(fā)生的情況的視覺和數(shù)值反饋，允許他們更好地理解和預(yù)測預(yù)期部件設(shè)計的行為，以便可以基本上消除制造的傳統(tǒng)、高成本的反復(fù)試驗的方法。計算機輔助工程模擬的使用便于在設(shè)計階段期間優(yōu)化部件設(shè)計、模具設(shè)計以及制造處理參數(shù)，在該設(shè)計階段中能夠以最低的成本和對進度的最小影響來容易地實施必要的改變。
CAE模擬技術(shù)在用于纖維增強部件的工程過程中的應(yīng)用包括(i)“注塑成型”制造工藝的模擬，其包括流體流和熱傳遞的計算；和(ii)壓力和強度(還可能有耐用性)計算，這些都是在宏觀水平上針對這些部件執(zhí)行的，以確定它們在外部負載下的功能機械性質(zhì)。兩種模擬類型都需要適當(dāng)?shù)牟牧夏Ｐ?，該材料模型描述包含浸入的纖維的聚合物材料在液態(tài)以及固態(tài)下的特性。
宏觀水平上的長度尺度由部件的幾何形狀的線性尺寸(整體大小、壁厚度等)確定，該線性尺寸通常在幾mm到cm的范圍內(nèi)變化。計算單元的尺寸必須以足夠的精度分解宏觀長度尺度，因此它們通常比最小宏觀尺寸小到一個量級。由于浸沒在短纖維增強部件中的纖維的典型尺寸低于宏觀計算單元的典型尺寸一個或者兩個量級，因此通過統(tǒng)計方法來描述與宏觀材料行為的建模有關(guān)的纖維性質(zhì)。對于短纖維增強材料，相關(guān)的宏觀性質(zhì)是(a)體積濃度，其關(guān)于整個零件通常(近似)恒定，和(b)在每個計算單元內(nèi)的纖維取向(FO)的局部分布，其通?？绮考缀涡螤蠲黠@變化。(在詳細描述的1.1和1.2節(jié)討論此問題的進一步細節(jié))。
通過相應(yīng)分布函數(shù)的低階(即二階和四階)矩提供局部FO的統(tǒng)計分布的簡化的、出于實際目的的適當(dāng)描述。由于它們的數(shù)學(xué)結(jié)構(gòu)，將這些矩表示為(分別為二階和四階的)取向張量。在用于纖維增強部件的CAE模擬的構(gòu)架內(nèi)，需要四階張量來預(yù)測纖維增強材料在宏觀水平上的流變和機械性質(zhì)，因為這些是四階張量性質(zhì)。二階FO張量是具有單位跡的實值對稱3×3矩陣，因此其9個分量中只有5個是獨立的。通過(全)對稱將四階FO張量的獨立分量的數(shù)量從34＝81減小到15。
以其矩的形式描述FO分布的數(shù)學(xué)模型通過根據(jù)閉合關(guān)系以二階張量的形式近似計算四階取向張量而被顯著地簡化。閉合關(guān)系以函數(shù)的形式提供這種計算方案的數(shù)學(xué)描述，并且如果閉合關(guān)系僅僅在特定假設(shè)下近似有效，則將相關(guān)計算過程稱為“閉合近似”。只使用二階FO張量和某些閉合近似的方法導(dǎo)致“Folgar-Tucker”型的模型以模擬在模具填充過程期間二階FO張量在時間和空間內(nèi)的演化(詳見詳細描述的2.2到2.5節(jié))。
文獻“Glass fibre orientation within injection moulded automotivepedal Simulation and experiment studies，B.R.Whiteside等人，Plastics，Rubber and Composites，2000，卷29，No.1”公開了用于對增強熱塑料物件內(nèi)的纖維取向分布進行建模的方法，其用于非對稱熱塑料流，并且分析包含纖維取向預(yù)測算法。該軟件使用由線性三角單元組成的二維有限元網(wǎng)格來近似三維模型。使用廣義Hele-Shaw近似和Folgar-Tucker方程的變型來計算流場以計算纖維取向。在穿過每個單元的“厚度”的19個夾層(laminate)上使用有限差分技術(shù)來進行纖維取向、溫度和粘度計算以產(chǎn)生三維解。然而，重要的是應(yīng)該注意到由于該模型不能模擬在平面之外方向上的速度分量(Hele-Shaw近似的限制)，因此不能將該系統(tǒng)描述為真實的三維。當(dāng)適用于真實三維模擬時，在該文獻中描述的方法將導(dǎo)致不穩(wěn)定以及過多處理能力消耗的模擬，該模擬不能在例如開發(fā)工程師的工作場所中使用。
發(fā)明內(nèi)容
在此背景下，本發(fā)明的目的是提供用于在過程模擬中確定非球形粒子的取向分布的方法，在該過程中，用包含大量非球形粒子的懸浮液填充模具型腔，與現(xiàn)有技術(shù)方法相比，該方法更加穩(wěn)定并且計算強度較小。
根據(jù)權(quán)利要求1，通過提供用計算機實現(xiàn)的方法來實現(xiàn)該目的，該方法用于通過使用模擬注塑成型過程的模擬模型來計算非球形粒子在宏觀水平上的取向統(tǒng)計，在該過程中，形成模擬區(qū)域的幾何形狀的至少部分的模具型腔被由包含大量非球形粒子的溶劑形成的懸浮液填充或部分地填充，其中提供了模擬區(qū)域的幾何形狀的數(shù)字表示或者計算機模型，并且其中通過對模擬區(qū)域的至少部分進行細分或者離散化來形成具有多個計算單元的網(wǎng)格，該方法包括步驟a)指定邊界條件；b)設(shè)置初始條件；c)針對模擬區(qū)域的至少部分單元求解質(zhì)量、動量和能量的平衡方程以獲得宏觀水平上的流體流、熱流和質(zhì)量傳遞；和d)至少部分基于所解的平衡方程的結(jié)果來求解非球形粒子取向動力學(xué)方程，以由此確定宏觀水平上非球形粒子取向的變化作為空間和時間的函數(shù)。這里，對于步驟d)，可能僅僅針對包含懸浮液的計算單元求解粒子取向方程。
優(yōu)選地，步驟c)還包括(cc)至少部分基于所解的平衡方程的結(jié)果確定流體或者懸浮液的更新的自由表面或者流前沿(flow front)，其中該自由表面將型腔中由懸浮液填充的單元與空單元分開。優(yōu)選地，步驟(cc)還包括根據(jù)所更新的流前沿更新邊界條件。同樣優(yōu)選地，本發(fā)明的方法還包括步驟e)通過確定模具型腔是否被懸浮液的模擬注射填滿來確定是否完成了模擬注塑成型過程；和f)重復(fù)步驟c)、cc)和d)，直到完成模擬注塑成型過程為止。
根據(jù)權(quán)利要求6，通過提供用于在過程模擬中描述非球形粒子的統(tǒng)計分布取向的方法，也實現(xiàn)了本發(fā)明的目的，在該過程模擬中，用由包含大量非球形粒子的溶劑形成的懸浮液填充模具型腔，該方法包括(1)提供定義型腔的幾何形狀的三維計算機模型；(2)指定邊界條件；(3)基于該模型對解域進行離散化以形成具有多個單元的網(wǎng)格；(4)針對解域的至少部分，求解能量和流方程；(5)計算各自單元中的流和溫度條件作為時間的函數(shù)；(6)計算非球形粒子取向的變化；和(7)將在各自單元中的非球形粒子取向的統(tǒng)計分布描述為時間的函數(shù)。
本發(fā)明的方法以顯著降低的計算強度利用懸浮液流來預(yù)測遍及部件的纖維取向分布和熱機械性質(zhì)。開發(fā)工程師可以使用模擬的結(jié)果，并能夠由此優(yōu)化與纖維取向相關(guān)的產(chǎn)品并且因此改進該物件的強度和形狀穩(wěn)定性。
根據(jù)詳細描述，根據(jù)本發(fā)明的方法和裝置的進一步的目標(biāo)、特征、優(yōu)勢和性質(zhì)將變得明顯，該方法和裝置用于確定在過程模擬中的非球形粒子的取向分布，其中在該過程中，用包含大量非球形粒子的懸浮液填充模具型腔。
在以下本說明書的詳細部分中，將參考附圖所示的示例性實施例更詳細地說明本發(fā)明，其中 圖1是穿過包括模具的注塑成型機的圖形表示的橫截面視圖； 圖2是概括根據(jù)本發(fā)明的實施例的注塑成型過程的基本處理步驟的上層流程圖； 圖3是進一步詳細圖解圖2的流程圖的步驟5的流程圖； 圖4是纖維增強塑料物件的顯微圖； 圖5圖解了在本發(fā)明的實施例中使用的纖維的模型； 圖6是圖解特征多項式的二次和三次以及線性和常數(shù)項的曲線圖； 圖7是圖解矩陣的三個特征值的圖解； 圖8示出了對應(yīng)于“極大對稱”情況的特殊示例； 圖9示出了四面體形體積的扭曲形式； 圖10示出了根據(jù)本發(fā)明的實施例的FO矩陣(即2階FO張量)的相空間集合MFT的總體結(jié)構(gòu)的圖； 圖11是圖解根據(jù)本發(fā)明的實施例的算子分裂過程是以簡化形式的流程圖； 圖12是圖解時間步長方法的細節(jié)的流程圖；和 圖13是圖解根據(jù)本發(fā)明的實施例的使用跡尺度改變的相空間投影過程的流程圖。
具體實施例方式 圖1圖示性地示出了注塑成型機1。該注塑成型機設(shè)置有螺桿2，其被饋送設(shè)置在漏斗3中的聚合物粒料。聚合物粒料通過螺桿2和加熱單元4的作用而轉(zhuǎn)變成粘滯團，該粘滯團在高壓下被推動到模具6中的模具型腔5中。并且成型機和注塑成型制造周期在本領(lǐng)域內(nèi)是眾所周知的，在此不詳細說明。使用成型機1，可以制造非纖維增強和纖維增強塑料部件。
可以根據(jù)圖2所圖解的過程在計算機上進行注塑成型過程的數(shù)值模擬。
一般所認(rèn)為的模擬的主要步驟如下 -步驟1，提供模擬區(qū)域的幾何形狀的數(shù)字表示； -步驟2，網(wǎng)格劃分(enmeshment)，其將計算區(qū)域細分成很多小單元，這些小單元是用于將微分方程離散化(利用不同求解算法)的基礎(chǔ)，并且以此方式找出要模擬的物理現(xiàn)象的解； -步驟3，將對于不同的材料域的必要物理數(shù)據(jù)附加到模擬模型中(數(shù)據(jù)庫或者資料庫)； -步驟4，指定針對模擬項目的邊界條件； -步驟5，模擬注塑成型過程(以下將更詳細地聲明該步驟)；和 -步驟6，在計算機(諸如工作站)的顯示器上將結(jié)果顯示為圖形或者數(shù)值表示。
在圖3的流程圖中圖解了當(dāng)模擬纖維增強物件的注塑成型過程時的步驟5的細節(jié)。在過程的這個部分，使用數(shù)值算法來解熱流、流體流以及壓力和張力的微分方程 -步驟1，設(shè)置熱物理材料性質(zhì)和流前沿的初始條件； -步驟2，使用質(zhì)量、能量和動量方程的守恒來求解整個域的熱方程和所有流體單元的流方程； -步驟3，在此步驟中移動流前沿并且根據(jù)新的流前沿采用邊界條件； -在步驟4中，根據(jù)在之前的步驟中獲得的流速計算纖維取向和輸運。
以下將更詳細地說明該步驟。將初始條件應(yīng)用于新填充的單元。僅考慮包含流體材料的單元； -在步驟5中，計算例如化學(xué)反應(yīng)的附加量，并且驗證單元是否凝固； -在步驟6中，驗證模具中的注塑成型過程是否結(jié)束；如果沒有結(jié)束，則模擬繼續(xù)下一時間步長并且該過程返回到步驟2； -在步驟7中，計算在從模具脫模之后的物件的性質(zhì)。
使用從注塑成型模擬獲得的信息計算在模具外的熱物理材料溫度，即溫度、收縮、翹曲等。
圖4是纖維增強塑料物件的顯微圖，其中可以看到在注塑成型過程結(jié)束之后的纖維取向。最終產(chǎn)品中的纖維取向很大程度上取決于注塑成型過程期間的熱塑性團的流型。不是針對每個單個纖維精確地確定纖維取向，而是通過分布函數(shù)進行描述。
在進入關(guān)于纖維取向計算的細節(jié)之前，提供短纖維熱塑性熔體的基本流體機械方面的概述。
1.短纖維熱塑性塑料熔體的基本流體機械方面 聚合物團包括分散于其中的大量的短纖維，使得所制造的部件將由短纖維增強熱塑性塑料制成。在熔化狀態(tài)中，即當(dāng)溫度足夠高使得熱塑性塑料基體是液體時，塑料熔體和分散的纖維的混合物組成復(fù)雜流體，該復(fù)雜流體在流體動力學(xué)和物理學(xué)的術(shù)語中通常被稱為粒子懸浮液。一般而言，這種懸浮液由兩種不同的相組成(i)溶劑，在我們的情況中對應(yīng)于沒有分散纖維的熔化的塑料材料，以及(ii)粒子相，其由所有浸沒在該溶劑中的纖維組成。
1.1纖維幾何形狀和材料性質(zhì) 在下文中將擁有旋轉(zhuǎn)對稱軸的球形粒子(如圖5所示)稱為纖維，并且除了另外明確指出以外，同義地使用術(shù)語粒子和纖維。通過纖維的長度l和直徑d以及從這兩個量得到的縱橫比ra＝l/d來表征纖維的幾何形狀。分散在短纖維增強熱塑性塑料材料中的纖維(例如，碳或者玻璃類型)的這些量的典型值例如是l≈0.5mm、d≈5μm以及ra≈100。通常這些值以相應(yīng)的大縱橫比值在l＝0.1到1mm、d＝1到10μm的范圍內(nèi)變化。
對于浸沒在纖維增強塑料中的大部分類型的纖維來說，纖維材料的質(zhì)量密度可與懸浮液塑料熔體的密度相比。塑料熔體本身在模具填充階段期間的典型溫度下具有相當(dāng)高的粘性。通過結(jié)合這兩方面，當(dāng)描述在根據(jù)優(yōu)選實施例的方法中的塑料熔體內(nèi)的纖維運動時忽略浮力效應(yīng)和慣性效應(yīng)。
關(guān)于纖維濃度，假定貫穿懸浮液的纖維的空間均勻分布，并且因此認(rèn)為濃度不變?？梢愿淖兊男Ч梢允?i)已經(jīng)在澆口處存在的不均勻分布的純粹對流輸運，或者(ii)在模具填充階段期間的剪切誘導(dǎo)粒子遷移效應(yīng)的出現(xiàn)。
如果(i)或者(ii)組成顯著效果，則用于描述浸沒在熱塑性塑料熔體中的短纖維的懸浮液流的模型有必要地需要是“兩相流”模型，該模型將粒子相和溶劑相的流彼此耦合并且允許不均勻粒子濃度。然而，情況并非如此?？梢詫F(xiàn)象(ii)理解為一種擴散過程，其中擴散常數(shù)與相對粒子大小的平方成比例，相對粒子大小定義為懸浮粒子的實際大小(即球形粒子情況中的纖維長度或者半徑)與流型腔的典型尺寸(例如，壁的厚度)的比。由于即便對于通過注塑成型由纖維增強塑料制造的薄壁部件而言，相對粒子大小通常為0.1到0.01的量級，因此可以完全合理地認(rèn)為在這種情況中(ii)完全可以忽略。關(guān)于(i)，可以想象由于在螺桿中的熔體的準(zhǔn)備過程，在澆口處可能存在不均勻濃度輪廓。然而，針對典型短纖維增強塑料部件的實際實驗性觀察的纖維濃度在部件內(nèi)的幾乎各個地方僅顯示非常小地偏離常數(shù)值，這可能是均勻纖維分布的假設(shè)的最佳論據(jù)[20]。
1.2纖維對流變性質(zhì)的影響 熱塑性塑料材料顯示非牛頓流行為。在沒有纖維的純塑料熔體的情況中，可以通過標(biāo)量粘度函數(shù)來對材料的流變性質(zhì)進行建模，標(biāo)量粘度函數(shù)取決于溫度以及流(廣義牛頓流體)的狀態(tài)變量。盡管已知它們僅覆蓋非牛頓流性質(zhì)(如例如剪切稀釋(shear thinning))的有限范圍，但是已經(jīng)證明它們對于注塑成型模擬的目的是有用的并且令人滿意地精確。此類模型也用于根據(jù)本發(fā)明的優(yōu)選實施例的方法中。
向熱塑性塑料基體材料加入纖維產(chǎn)生的機械性質(zhì)在固態(tài)下是強烈各向異性的，并且該機械性質(zhì)很大程度上取決于嵌入的纖維所指向的方向的局部分布。原則上，如果材料在液(即熔化)態(tài)下，則也存在各向異性材料行為。為了解決這種各向異性，不得不用粘度張量取代上述標(biāo)量粘度。已經(jīng)存在若干研究對纖維懸浮液的材料行為的兩種建模方式進行比較。對于如在模具填充中遇到的大部分流情形，已經(jīng)發(fā)現(xiàn)通過兩種類型的模型(即標(biāo)量和張量模型)所預(yù)測的充填模式顯示出幾乎沒有差別(參見例如[21]、[22])。因此當(dāng)熱塑性塑料熔體包含纖維時，忽略各向異性粘度效應(yīng)并且在這里也使用簡單的廣義牛頓模型。這等價于忽略分散的纖維的取向?qū)θ垠w的流變性質(zhì)的影響。粘度取決于纖維濃度，但是由于假設(shè)纖維濃度不變(見上)，因此這個方面僅作為先驗已知的材料參數(shù)輸入，該參數(shù)有助于熱塑性塑料熔體的有效材料性質(zhì)。
由于這種途徑，根據(jù)優(yōu)選實施例的方法使用流和纖維取向的計算的部分解耦盡管流速局部地影響分散纖維的取向，但是纖維取向?qū)α鞯挠绊懣珊雎?。因此獨立于纖維取向的計算進行流的計算。以此方式，熔體的局部流速作為一組外部系數(shù)進入用于纖維取向計算的模型。
2.Folgar-Tucker模型 2.1 Jeffry方程 盡管懸浮在聚合物熔體中的纖維在它們的縱橫比ra＝l/d較大的意義上說是細長粒子，但是它們足夠短以使得在溶劑的局部流場中作用在其上的機械力不能夠引起任何實質(zhì)的變形。因此這里將單個纖維建模為細長、旋轉(zhuǎn)對稱的剛體，其取向由沿著旋轉(zhuǎn)對稱軸指向的單位向量p給出。兩個向量±p都表示粒子的相同取向狀態(tài)。
圖5圖解了旋轉(zhuǎn)橢圓形的剛性粒子，該粒子的運動受粒子附近的速度向量場U影響。除了表征粒子的幾何形狀的量(長度l和直徑d)之外，示出了取向向量p。由相應(yīng)箭頭的方向和長度指示流速U的方向和大小。盡管速度向量的方向都一樣，但是其長度不同，指示粒子附近的流速不是恒定的。然而，由于從底部到頂部向量的長度有恒定的增加，因此流速變化的局部量(即速度梯度)看起來相同。在粒子表面的每個點處，認(rèn)為粒子精確地以局部流速(“無滑動(no-slip)”邊界條件)移動。如果粒子周圍的流速是恒定的，則這將導(dǎo)致簡單的平移運動，即純粹對流輸運。否則，在存在速度梯度的情況下，如虛線箭頭所指示的，粒子還進行旋轉(zhuǎn)運動。總而言之，剛性粒子進行的最普通類型的運動是以粒子周圍的“平均速度”進行的平移，結(jié)合圍繞其質(zhì)心的旋轉(zhuǎn)，該旋轉(zhuǎn)由描述粒子附近的流速與其平均值的局部偏差的速度梯度驅(qū)動。
已經(jīng)由Jeffery將以上所給出的定性描述實現(xiàn)為數(shù)學(xué)模型[2]，其假定粘性力是主導(dǎo)的，慣性力因此可忽略，并且在單個纖維運動期間由其掃過的區(qū)域上的局部流的變化較小。如上所說明的，所有這些假設(shè)都在分散在熱塑性塑料熔體中的短纖維的情況中得到了滿足。所假設(shè)的纖維周圍的流速變化較小意味著局部速度梯度張量
足以精確地描述該變化。該張量是3×3矩陣，根據(jù)流速向量的分量Uj關(guān)于空間中的點r＝(x1，x2，x3)T的笛卡爾坐標(biāo)xi的偏導(dǎo)數(shù)(即)計算該矩陣的元素。因此如 果r0是粒子質(zhì)心的空間坐標(biāo)，則通過一階Taylor展開式來具有可忽略的誤差地較好地近似在粒子附近的點r處和時間t處的流速U的值，并且可以把速度梯度看作局部恒定量，即可以假設(shè)適用于r0附近的所有點r。
在這種情況下，使用以下Jeffery方程計算由作為空間和時間坐標(biāo)的函數(shù)的其取向向量p(r，t)給出的纖維的瞬時取向狀態(tài) 這里以緊湊歐拉形式寫出該方程。方程(1)的左側(cè)(l.h.s.)的對流(或者材料)導(dǎo)數(shù)描述了純粹纖維取向FO、流的局部速度場U(r，t)中的輸運，而右側(cè)(r.h.s.)對由有效局部速度梯度張量
所驅(qū)動的粒子旋轉(zhuǎn)運動進行建模。在無限細長纖維(即縱橫比ra→∞)的理想特殊情況下，張量
等同于速度梯度張量
在纖維具有有限縱橫比(0＜ra＜∞)的一般情況下，由以下項定義張量
該項包含纖維幾何形狀參數(shù)λ＝(ra2-1)/(ra2+1)該參數(shù)對在旋轉(zhuǎn)對稱橢圓粒子的情況下粒子幾何形狀如何影響纖維的旋轉(zhuǎn)運動進行編碼。
利用唯一定義的分解將速度梯度張量分解成由剪切速率和渦度張量給出的其對稱和不對稱部分，以該可選方式將Jeffery方程寫為 (“剪切速率”)， (“渦度”)。
剪切速率和渦度張量經(jīng)由以下恒等式關(guān)聯(lián)到有效速度梯度張量 將這些恒等式代入(1)并且考慮產(chǎn)生以下的可選形式Jeffery方程 仍然可以通過使用恒等式來重寫(1d)的右側(cè)的第一項來獲得Jeffery方程的另一形式，該恒等式將渦度張量關(guān)聯(lián)到流速向量場的旋度rotU(見例如[3])。
在粒子的幾何形狀是旋轉(zhuǎn)對稱但不是橢圓的情況下，方程(1)的形式以及有效速度梯度的定義保持不變，但是必須改變幾何形狀參數(shù)λ的公式，使得用有效縱橫比取代如上定義的“幾何形狀”縱橫比ra＝l/d，需要根據(jù)實驗確定該有效縱橫比的合適值。在這種意義上，可以將幾何形狀參數(shù)λ看作材料參數(shù)。在非對稱非球形粒子的情況下，需要用包含三個幾何形狀參數(shù)的三個Jeffery型方程的耦合系統(tǒng)取代方程(1)[3]。
2.2纖維取向的宏觀分布 在宏觀水平上執(zhí)行針對短纖維增強塑料的注塑成型模擬的情況下，因為單個纖維的尺寸小，因此包含在單個計算單元(即計算區(qū)域的體積單元)中的纖維的數(shù)量非常大。這一點由圖4所示的顯微解，該顯微圖表示從由短纖維增強塑料制成的部件得到的典型樣本。因此通過宏觀模型(即FO向量p的分布函數(shù)Ψ(p，r，t)(見[5]))描述包含在計算單元中的材料的局部FO狀態(tài)。我們注意到盡管假設(shè)纖維濃度是均勻的，但是分布函數(shù)ψ仍然經(jīng)由局部流速場U(r，t)及其梯度
參數(shù)地依賴于(t，t)坐標(biāo)。作為附加點，需要將說明纖維的相互作用、影響它們的取向的項包括在宏觀水平上的模型中。
經(jīng)由以下相應(yīng)的Fokker-Planck方程實現(xiàn)將作為單獨的、無相互作用的纖維的微觀模型的Jeffery方程轉(zhuǎn)換成產(chǎn)生許多相互作用纖維的FO統(tǒng)計的宏觀模型 其中以FO分布函數(shù)Ψ作為因變量。這里
作為Jeffery方程(1)的右側(cè)的簡寫，并且
和
是在三維歐幾里得(Euclidian)空間
的單位球體S2上定義的梯度和散度算子的符號。根據(jù)Folgar和Tucker的方法[4]，與局部速度梯度的有效剪切速率成比例的擴散系數(shù)Dr＝CIγeff的引入產(chǎn)生在濃縮懸浮液中纖維纖維相互作用的簡單模型。這里Gij是如(1b)中定義的剪切速率張量
的分量，并且平方根號下的項(使用愛因斯坦求和約定)是其自收縮的簡寫。無量綱、非負相互作用系數(shù)CI是懸浮液的材料參數(shù)。通常對于短纖維增強塑料，該材料系數(shù)具有在10-3到10-2范圍內(nèi)的小(正)值。我們注意到在不可壓縮的(以及幾乎不可壓縮的)流中，“隨機”擴散項(與表示“確定性”Jeffery動力學(xué)的項相比)的相對弱點可能是穩(wěn)定性問題的來源。
2.3纖維取向張量和Folgar-Tucker方程 通過Fokker-Planck方程(2)計算局部FO分布需要在針對流模擬區(qū)域的每個計算單元的單位球體S2上定義的PDE的數(shù)值解，這是對于“工業(yè)尺度”3D問題來說負擔(dān)不起的昂貴的任務(wù)。因此，Advani和Tucker[6]提出使用纖維取向張量，其被定義為分布函數(shù)的矩，并且因此用FO張量的一族矩方程取代Fokker-Planck方程。由于FO分布關(guān)于變量p的反轉(zhuǎn)對稱Ψ(-p，...)＝Ψ(p，...)(這反映±p的方向?qū)?yīng)于相同的取向狀態(tài)的事實)，因此所有奇數(shù)階的矩一致地消為零，以使得Ψ的矩展開式僅包含偶數(shù)階的元素。因此這個展開式的第一非平凡矩是第二個，其由下式給出
(以索引表示
符號
表示在球體S2的表面上進行積分，其中dS(p)是積分度量。第二矩
被稱為二階FO張量(或者FO矩陣)并且根據(jù)定義是實對稱3×3矩陣。由于根據(jù)
標(biāo)準(zhǔn)化FO分布(也根據(jù)定義)，因此當(dāng)滿足時，F(xiàn)O矩陣
具有明顯的性質(zhì)，即其跡等于1。展開式體系中的下一個非平凡矩是四階FO張量
其定義為
(以索引表示
張量
是完全對稱的并且額外地擁有各種標(biāo)準(zhǔn)化性質(zhì)由于p2＝1，因此相同索引的任意對的和總是產(chǎn)生
的相應(yīng)元素(例如，)，并且兩個相同索引對的和總是等于1(例如，)，因此
包含關(guān)于
的完整信息。
可以通過交換微分和積分運算來獲得針對每階的矩的上述矩方程族，即
用Fokker-Planck方程(2)的右側(cè)的項取代
并且解析地估算相應(yīng)的積分。如果將這個過程應(yīng)用于FO矩陣
，則獲得所謂的Folgar-Tucker方程(FTE)作為分級組織的一系列矩方程中的第一非平凡方程 如在Jeffery方程中，(5)的左側(cè)的對流導(dǎo)數(shù)描述由歐拉參考系中纖維的平移運動造成的局部FO狀態(tài)(由矩展開式的該階的FO矩陣
表示)的純粹輸運，而(5)的右側(cè)的前兩項對由局部有效速度梯度
所驅(qū)動的它們的旋轉(zhuǎn)運動進行建模。從Fokker-Planck方程(2)中的擴散項的存在產(chǎn)生(5)的右側(cè)的第三項。在FTE的水平上，其產(chǎn)生阻尼效果，該效果將FO矩陣拖向各向同性狀態(tài) 在其歐拉形式(5)中，F(xiàn)TE是對流反應(yīng)(convection-reaction)類型的一階PDE的耦合系統(tǒng)。從拉格朗日觀點來看，(5)是針對
的分量的ODE的耦合系統(tǒng)。
2.4FO矩陣的基本性質(zhì) 作為實對稱3×3矩陣，F(xiàn)O矩陣
擁有3個實特征值μk，μk的相應(yīng)特征向量為Ek，Ek形成
的標(biāo)準(zhǔn)正交基，即且Ek·El＝δkl，其中δkl是克羅內(nèi)克(Kronecker)符號，對于k＝l，該符號等于1，否則等于零。顯而易見的是，通過構(gòu)造，F(xiàn)O矩陣
的特征值μk是非負的并且(當(dāng)滿足p2＝1時)它們的和(其與
的跡相同，見上)總是等于1(即)。
對于任意單位向量p0，特殊FO矩陣表示所謂的單軸取向狀態(tài)，其中100％的纖維在±p0給出的方向上取向。因此，p0的符號不參與定義這個特殊FO矩陣的雙積乘積。向量p0是
的特征值為1的特征向量，并且其相應(yīng)的特征向量位于與p0正交的平面中的剩余的兩個特征值為零，否則為任意。
可以其譜表示的方式來寫FO矩陣，該譜表示具有由FO矩陣的特征向量形成的單軸取向狀態(tài)
的加權(quán)和的形式，權(quán)重由特征值給出。這允許將特征值μk解釋為沿著相應(yīng)的特征向量Ek的方向取向的纖維的局部部分。在這個意義上，F(xiàn)O矩陣的譜數(shù)據(jù){μk，Ek}k＝1，2，3表示在懸浮液的小體積內(nèi)的纖維的局部宏觀取向狀態(tài)。
這作為FO矩陣和FTE的相空間的以下(數(shù)學(xué)形式)定義的動機[12]定義當(dāng)且僅當(dāng)實對稱3×3矩陣
是半正定并且其跡等于1時，該矩陣是FO矩陣。FTE的相空間MFT是所有FO矩陣的集合。
可以示出的是，集合MFT等于使用適當(dāng)?shù)貥?biāo)準(zhǔn)化的(但是否則為任意的)分布函數(shù)從類型(3)的矩積分得到的所有實對稱3×3矩陣的集合。最近由作者之一利用相空間MFT的拓撲和幾何性質(zhì)給出了其數(shù)學(xué)特征化(見[12])。由于FTE的數(shù)值積分結(jié)果的可解釋性需要因變量
在FO計算的所有步驟(即，也包括中間步驟)過程中具有如上所述的特殊譜性質(zhì)，因此確切地知道矩陣是否屬于集合MFT非常具有實際重要性。
2.5閉合問題 所謂的閉合問題源自這樣的事實，即在矩展開式的每一階，矩
的DE包含作為變量的下一個更高階的矩
。盡管可以通過簡單代數(shù)恒等式(即對一對兩個相等索引求和，見上)來根據(jù)
表示
但是反過來是不可能的，因此需要將
視為未知。在FTE中，由(5)的右側(cè)的
的存在而顯示出了閉合問題，這妨礙了系統(tǒng)的可解性，除非通過將
通過閉合+近似表示為
的函數(shù)來閉合該系統(tǒng)。將閉合近似應(yīng)用于FTE意味著由FO矩陣的某個合適的(一般為非線性)張量函數(shù)
來取代(5)的右側(cè)的精確(但未知的)四階FO張量
公知的示例是由以下公式定義的混合閉合[7] 盡管有一些眾所周知的缺點，但是由于該公式的代數(shù)簡單性以及數(shù)值魯棒性因而其是可接受的選擇[7]。將閉合
定義為在兩個較簡單類型的閉合之間的(凸面)插值二次閉合定義為 (即以索引表示為) 其在單軸取向分布的特殊情況中產(chǎn)生精確的結(jié)果，并且由下式給出線性閉合 (6b) 其在各向同性的情況下是精確的。由標(biāo)量取向因子提供在這些極端情況之間的插值權(quán)重 由于行列式det
是FO矩陣的不變量，因此對于標(biāo)量取向因子同樣成立。
我們通過首字母縮寫FTE-hyb來表示結(jié)合混合閉合近似(6)的FTE(5)的特殊變量。該變量非常具有實際意義，本發(fā)明的優(yōu)選實施例也是基于Folgar-Tucker模型的FTE-hyb變量。
2.6FTE作為微分代數(shù)系統(tǒng) 針對任意實對稱矩陣而形式上良好地定義(5)的右側(cè)。對于(5)的右側(cè)的性質(zhì)是否可以自動將解軌跡限制到域MFT這個問題，可以至少針對其精確形式(即無閉合近似)下的FTE來肯定地回答如果將精確四階FO張量
代入(5)的右側(cè)中，并且然后作為(5)的解計算
則結(jié)果與FO分布的精確二階矩相同，其通過使用具有相同的參數(shù)
和Dr的Fokker-Planck方程(2)的解Ψ估算動量積分(3)來獲得。通過這個論據(jù)，可以推斷，使用精確四階FO張量
可以將(5)的解寫為矩積分(3)，并且因此必然地將其軌跡限制到域MFT。
然而，如果應(yīng)用了閉合近似，則這個論據(jù)不再有效，如果必須在對完全分布函數(shù)沒有任何先前知識的情況下數(shù)值求解FTE，則應(yīng)用閉合近似總是必要的。以此方式，將FTE的解限制到相空間MFT的問題總是出現(xiàn)在所有的具有實際意義的問題中。
將FTE的解軌跡限制到域MFT的必要性給因變量
加入附加限制，可以依據(jù)其不變量的代數(shù)不等式來用公式表示該限制(見[12])。這些關(guān)于因變量的限制將FTE變成微分代數(shù)系統(tǒng)(DAS)并且需要在用于數(shù)值積分的程序中注意該限制。
3.FTE的相空間的數(shù)學(xué)特征化 可以根據(jù)FO矩陣的特殊性質(zhì)直接推導(dǎo)出相空間MFT的基本拓撲性質(zhì)，其中FO矩陣的特殊性質(zhì)概括如下 定理1相空間MFT是被限制到由跡條件定義的五維超平面的實對稱3×3矩陣的向量空間的有界凸子集。
MFT的凸度考慮到向MFT的投影映射的定義，其可以與任何適當(dāng)?shù)腛DE積分器組合以構(gòu)成用于FTE的合適的積分方法(見[9]的第IV.4章)?？梢酝ㄟ^對實對稱3×3矩陣的以下特征多項式的分析獲得FO矩陣的不變量代數(shù)特征化 系數(shù)和是矩陣的不變量。(在文獻中，有時由I1＝Sa、I2＝Ka和I3＝Da表示這些不變量。)可以根據(jù)這些不變量用公式表示FO矩陣的代數(shù)特征化，其根據(jù) 定理2當(dāng)且僅當(dāng)實對稱3×3矩陣
的跡Sa等于1并且其不變量Ka和Da是非負時，該矩陣是FO矩陣。
圖6通過對特征多項式Pa(μ)的二次和三次項以及線性和常數(shù)項的單獨分析而圖解了此表述正跡Sa明顯地提供正特征值μk的存在，而兩個不變量Ka和Da的非負值阻止了負特征值的存在。此外，由于矩陣
總是具有三個實特征值，因此條件Sa＞0，Ka≥0和Da≥0對于
的所有特征值為非負不僅是必要的而且是充分的。結(jié)合跡條件Sa＝1，這完成了以上定理的證明。
可以通過分別查看FO矩陣的對角和非對角元素來獲得相空間MFT的幾何形狀的概念，該相空間根據(jù)定理1是5D對象。首先我們注意到由
給出的對角元素總是非負的并且滿足跡條件∑kakk＝1。因此，不依賴于坐標(biāo)系的選擇，它們被限制到三角集合(見圖7)
這提供了必要條件，該必要條件將FO矩陣關(guān)于它們的對角元素特征化并且產(chǎn)生相空間集合MFT的“對角部分”的不變量表示。
可以通過為實對稱3×3矩陣的對角和非對角元素的三元組引入符號(x，y，z)和(u，v，w)，來獲得域MFT的“非對角”部分的正規(guī)描述，取任意(但固定)的對角三元組(x，y，z)∈Δa并且形式上定義屬于固定對角三元組的所有“容許的”非對角三元組的集合
如關(guān)于以上定理2所說明的，在代數(shù)上可以將集合N(x，y，z)特征化為同時滿足下面一對不等式的所有非對角三元組(u，v，w)的集合 在圖8中示出了對應(yīng)于“極大對稱”情況x＝y(tǒng)＝z＝1/3的特殊示例，其它情況(見圖9)對應(yīng)于這個四面體形的體積的扭曲形式。
為了補充不變量Ka和Da為非負值的限制(9)和(10)，我們注意到對于任何具有正跡Sa＞0的實對稱3×3矩陣，也根據(jù)以上由和來限制這些不變量，以便對于FO矩陣(Sa＝1)來說，總是由0≤Ka≤1/3和0≤Da≤1/27限制它們的值。
通過引入集合Δa和N(x，y，z)而獲得的(形式)纖維化有助于得到MFT的總體結(jié)構(gòu)的圖畫(見圖10)可以通過局部地將可容許的非對角三元組的集合N(x，y，z)附連到其“基點”(x，y，z)∈Δa的過程使單個“纖維”可視化，并且該基點貫穿整個三角集合Δa的隨后的變化覆蓋整個相空間。
在本節(jié)概括的結(jié)果清楚地顯示了相空間MFT是復(fù)雜的數(shù)學(xué)對象。根據(jù)定理2中敘述的嚴(yán)密數(shù)學(xué)結(jié)果，將被自然地定義為對稱3×3矩陣的實向量空間中的演化方程的FTE的任何解軌跡
限制到相空間域必然地意味著結(jié)合單位跡條件的不變量不等式(9)和(10)的滿足。這個事實不可避免地將FTE轉(zhuǎn)變成DAS。與單位跡條件不同，在存在任何目前公知的閉合近似的情況下，一般不保持不變量不等式(9)和(10)的有效性，而在不存在針對四階FO張量的任何閉合近似的情況下，單位跡條件的滿足被“構(gòu)建到”精確FTE中，并且在對于較大類的閉合近似的相當(dāng)普遍的先決條件下，單位跡條件仍然有效(見第4節(jié))。(同樣地，在“精確”FTE的情況下，可以寧可通過如2.6節(jié)給出的間接推論而不是通過FTE本身的代數(shù)結(jié)構(gòu)推導(dǎo)(9)和(10)的有效性。)在學(xué)術(shù)工作中經(jīng)常示出的簡單例子中通常忽略這些數(shù)學(xué)事實。
然而，任何數(shù)值求解FTE的過程有必要包含閉合近似作為其基本元素。因此，任何適合于處理更復(fù)雜的情形(如通常在工業(yè)應(yīng)用中出現(xiàn)的情形)的模擬過程必須通過適當(dāng)?shù)目刂七^程來對此進行解釋，該控制過程將解軌跡限制到其理論上可容許的域。因此對于嚴(yán)格的工業(yè)應(yīng)用來說，將FTE作為微分代數(shù)系統(tǒng)的適當(dāng)處理是強制性的。
在7.7節(jié)中，提出了提供適當(dāng)類型的不變量控制的過程。已經(jīng)在優(yōu)選實施例中實施了該過程。
4.跡守恒和穩(wěn)定性的問題 將閉合近似應(yīng)用于FTE意味著以FO矩陣的某(一般為非線性)張量函數(shù)
取代(5)的右側(cè)的精確(但是未知)四階FO張量
(即數(shù)學(xué)表示為進一步的細節(jié)見以上關(guān)于閉合問題的章節(jié))。取決于閉合的選擇，四階張量
具有精確張量
的較小或者較大量的性質(zhì)。
一個被所有合理的閉合近似滿足的關(guān)于
的對稱性質(zhì)的特殊要求是所謂的正交各向異性對稱性，其由這組等式給出 即，要求4階張量
關(guān)于第一和第二索引對(ij)和(kl)以及兩對的互換對稱。如果假設(shè)
是由具有正交各向異性對稱性質(zhì)(11)的閉合近似所增廣的(5)的解，則我們可以形式地推出其跡的以下微分方程(縮寫DE)(見[17]) 在(12)的推導(dǎo)中還不使用跡條件因為應(yīng)該通過對該方程的分析研究FTE的解的跡的穩(wěn)定性和守恒。如果將閉合近似應(yīng)用于FTE，則后者還滿足標(biāo)準(zhǔn)化條件 并且跡的DE(12)簡化為 由于擴散參數(shù)Dr根據(jù)定義是非負的，因此如果閉合近似滿足條件(11)和(13)，則可以推出(在Dr＞0的情況下)由跡條件定義的超平面是穩(wěn)定積分流形或者(在Dr＝0的情況下)跡是FTE的第一積分。(應(yīng)該注意的是，精確四階FO張量
根據(jù)定義是完全對稱的并且滿足并且精確FO矩陣
根據(jù)定義是對稱的并且滿足跡條件。) 以上的考慮顯示了取決于近似張量函數(shù)
和精確張量
共同具有的性質(zhì)的量，將閉合近似應(yīng)用于FTE如何影響跡條件(其為施加到FO矩陣的不變量的最簡單的條件)的有效性。
在混合閉合(6)的特殊情況下必須認(rèn)真考慮跡穩(wěn)定性，其中混合閉合僅擁有對稱性質(zhì)(11)但是不能滿足標(biāo)準(zhǔn)化條件(13)。在這種閉合近似的情況下，跡的DE采用以下形式(另見[17])
其具有前因子
該前因子偏離對于滿足(13)的所有閉合有效的較簡單的表達式6Dr，這反映了這樣的事實，即混合閉合未能滿足(13)。
在下文中示出了在某些條件下，在前因子
中出現(xiàn)的附加項可能使該項變?yōu)樨?。這意味著盡管超平面仍然是具有混合閉合的FTE的積分流形，但是其可能變得局部不穩(wěn)定。
可以通過以下論證的四階段過程來推出前因子
變?yōu)樨摰目赡苄? (i)在不可壓縮流場(即)的特殊情況下，前因子采用簡化形式
(ii)在細長纖維的情況下，可以總是假設(shè)λ＝1，并且如果由
描述的局部取向狀態(tài)是近似單軸(即其特征值μk之一接近1)，則由于
的行列式接近0，因此標(biāo)量取向因子的值也近似為1。以此方式，獲得簡化表達式
其在特定的假設(shè)下接近
的精確值。
(iii)通過代入Dr＝CIγeff并且將譜表示代入收縮獲得前因子的近似表達式的經(jīng)修改的形式
(iv)在最后一步中，對于有效剪切速率標(biāo)準(zhǔn)化剪切速率張量
的分量。標(biāo)準(zhǔn)化的剪切速率張量具有與
相同的特征向量、零跡并且因此(與
一樣)具有負特征值。由于根據(jù)構(gòu)造
的分量是1級的，因此張量
必然具有至少一個具有1級的絕對值的負特征值。
在到目前為止針對前因子
的近似表達式的推導(dǎo)中所假設(shè)的各種環(huán)境下，這個事實是
的期望的估計的關(guān)鍵。使用標(biāo)準(zhǔn)化的剪切速率張量
可以用以下形式寫出該近似表達式
考慮到在實際的模擬中，相互作用系數(shù)CI是小正數(shù)(通常為0＜CI＜10-2，另見1.2節(jié))，可以看出，如果以具有主特征值μj≈1的FO矩陣估算，
將變?yōu)樨摚撎卣髦稻哂邢鄳?yīng)的特征向量Ek，該特征向量沿著剪切速度張量的最大負特征值的特征方向指向，因為在這種情況下，可以估計并且因此得到
可以以嚴(yán)密的方式用公式表示以上論證。以此方式，可以給出數(shù)學(xué)證明，即針對任何不可壓縮流和相互作用系數(shù)的值0≤CI＜λ，具有混合閉合的FTE的解的跡在相空間MFT的某些區(qū)域中變得局部不穩(wěn)定，因為前因子
的表達式(15b)在這些區(qū)域中不可避免地變?yōu)樨摗?br>
因此，很清楚有必要在FTE的數(shù)值積分過程中注意跡穩(wěn)定性問題。
5.FTE的數(shù)值積分一般方面 用于FTE(包括任何類型的閉合近似)的數(shù)值積分的適合的方法的選擇取決于方程類型的數(shù)學(xué)分類以及與FTE的特定代數(shù)形式相關(guān)的方面。如在以下的段落中所詳細說明的，閉合FTE的一般結(jié)構(gòu)顯示出其構(gòu)成FO矩陣的作為因變量的元素aij(2)的對流反應(yīng)(convection-reaction)型的雙曲偏微分方程(PDE)的耦合系統(tǒng)。
從方程(5)中給出的FTE的形式開始，可以通過明確地代入(5)的左側(cè)的實質(zhì)導(dǎo)數(shù)(material derivative)和(5)的右側(cè)的象征閉合近似的數(shù)學(xué)符號來重寫這個方程，其產(chǎn)生閉合FTE的等價公式 FTE的這種形式已經(jīng)揭示了其數(shù)學(xué)結(jié)構(gòu)的大部分左側(cè)由輸運或者對流類型的簡單一階偏微分算子組成，其中局部流速U(r，t)控制作為方程(16)的因變量的FO矩陣的元素aij(2)的解耦、純粹對流輸運。此外，(16)的右側(cè)的各種項的代數(shù)結(jié)構(gòu)顯示，右側(cè)(正如FO矩陣本身，并且如由數(shù)學(xué)一致性所要求的)組成二階對稱張量函數(shù)，在以下方程(16)的簡要(但是等價)形式中將該函數(shù)表示為
盡管FO矩陣
的右側(cè)的顯式依賴性是線性的，但是由閉合近似的函數(shù)形式所確定的其隱式依賴性的性質(zhì)一般而言是非線性的。因此除非選擇了線性閉合近似，否則需要把函數(shù)
作為因變量
的非線性函數(shù)考慮。此外，可以認(rèn)識到的是，右側(cè)函數(shù)
不可避免地導(dǎo)致針對FO元素aij(2)的單獨方程的相互耦合，除非有效速度梯度
是對角矩陣并且閉合函數(shù)
具有非常特殊的代數(shù)結(jié)構(gòu)。
在通過方程(1a)解出了
對實際速度梯度張量
和纖維幾何形狀參數(shù)λ的依賴性并且以有效剪切速率的形式用定義表達式Dr＝CIγeff重新替代旋轉(zhuǎn)擴散參數(shù)之后，可以用以下形式重寫方程(17a) 閉合FTE的這種形式顯式地揭示了右側(cè)關(guān)于恒定模型參數(shù)λ和CI的依賴性，但是仍然在對FO矩陣的顯式和隱式依賴性之間進行區(qū)分，以及在對速度梯度張量
的顯式依賴性和由有效剪切速率公式引起的對
的隱式依賴性之間進行區(qū)分(詳見關(guān)于Fokker-Planck方程的章節(jié))。在文獻中，有時使用索引表示以分量形式寫方程(17b)，即 經(jīng)常使用精簡符號找到該方程的更加簡化的形式該精簡符號隱藏了所有的參數(shù)依賴性并且不在分量函數(shù)Fij(2)(...)對FO矩陣或速度梯度的顯式和隱式依賴性之間進行區(qū)分。從索引表示回到在方程(17a/b)中使用的完全張量符號，獲得了這些方程的簡化形式 6.FTE的數(shù)值積分“算子分裂” 存在各種用于對流反應(yīng)型的耦合PDE系統(tǒng)的數(shù)值積分方法。一種已經(jīng)被證明尤其在非線性右側(cè)函數(shù)F(...)的情況下魯棒并且靈活的方法從等價地將PDE重整為開始，并且通過將由線性微分算子
和函數(shù)F(...)確定的“算子”所給出的右側(cè)的兩項分開處理而進行。在數(shù)學(xué)文獻中，這種方法被稱為[23-37]“分步法(method of fractional steps)”、“分裂法(splitting method)”或者簡單地“算子分裂(operator splitting)”(關(guān)于可替選方法，參見[28-31])。通過“算子分裂”對具有形式的方程進行的數(shù)值積分利用兩個單獨方程和的數(shù)值(或者在某些情況下甚至是解析的)解，通過將右側(cè)的兩個算子之一設(shè)為零而從完整方程獲得所述兩個單獨方程。第一個方程等價于齊次對流方程第二個方程組成常微分方程(ODE)的系統(tǒng)，其一般而言是耦合并且非線性的。在優(yōu)選實施例中使用這種類型的方法。使用方程(17a/b)的符號，將分裂方法應(yīng)用于閉合FTE導(dǎo)致以下偏微分方程 (→在右側(cè)具有零矩陣
的“對流步驟”) 和 (→“旋轉(zhuǎn)步驟”) 將其視為在“算子分裂”方法的構(gòu)架內(nèi)的獨立的子問題。由方程(18a)建模的物理效果是由于根據(jù)流速的流體質(zhì)量的純粹對流輸運在模具型腔的被填充域內(nèi)的FO統(tǒng)計(如由FO矩陣的元素所編碼的)的空間再分布，而方程(18b)完全忽略對流輸運的影響并且單獨考慮由于由局部速度梯度驅(qū)動的纖維的旋轉(zhuǎn)運動以及纖維間的相互作用造成的FO分布的局部變化速率。在下文中，描述若干變型，其示出可以如何組合兩個子問題的(數(shù)值)解來產(chǎn)生完整方程(16)(或者其以簡要表示等價形式(17a/b))的近似解。
6.1“簡單算子分裂” 使用最簡單的方式(通常表示為“簡單算子分裂”)來獲得實際上根據(jù)兩個偏微分方程求解的方程的近似解，以連續(xù)順序求解第二個方程，將從第一個方程得到的中間解作為第二個方程的輸入(即初始值)。為了詳細說明，我們通過使用具有標(biāo)識和的模型方程來簡化符號。
流求解器(即對懸浮液的流進行建模的軟件)通過在位于包含在覆蓋模具型腔、澆注系統(tǒng)和澆口的總體計算區(qū)域的填充部分Ω(n)內(nèi)的空間中的離散點rm附近的所有計算單元中、在離散時間瞬時tn處求解針對質(zhì)量、動量和能量的離散輸運方程來產(chǎn)生流體的狀態(tài)變量，其中包括流速U(r，t)和其梯度張量
流求解器在時間t0＝0處開始，并且使用步長Δtn+1＝tn+1-tn從時間tn進行到tn+1。該流計算意味著在時間tn+1處的域Ω(n+1)的計算取決于其之前的狀態(tài)Ω(n)和新的速度場U(n+1)?？梢酝ㄟ^使用流體體積(VOF)方法來完成該計算。計算新的流前沿和新的域Ω(n+1)產(chǎn)生在時間tn+1處的Ω(n+1)的單元中的所有狀態(tài)變量的新值。在相應(yīng)的步驟中，從在根據(jù)之前的計算步驟已知的“舊”域Ω(n)中定義的值開始，應(yīng)該通過PDE的數(shù)值解來計算位于新填充的域Ω(n+1)的點rm附近的所有計算單元中的時間tn+1處的新值 這通過以下三步過程完成，該過程組成了“簡單算子分裂”的一種可能的變型 1.擴展步驟從域Ω(n)的單元中的開始，計算在新域Ω(n+1)的所有單元中的初始值ym(n)。
2.對流步驟從擴展步驟所提供的初始值ym(n)開始，使用域Ω(n+1)中的由流求解器所計算的流速U(rm，tn+1)，通過具有步長Δtn+1的對流方程的數(shù)值解計算中間值
3.旋轉(zhuǎn)步驟從之前的對流步驟所提供的初始值
開始，使用域Ω(n+1)中的由流求解器提供的速度梯度
通過具有步長Δtn+1的ODE系統(tǒng)的數(shù)值解計算最終值 圖11圖解了這種算子分裂變型。
上述的“簡單分裂”變型產(chǎn)生完整方程的近似解，其步長具有的一階時間離散化誤差(即O(Δtn+1))，并且其用于根據(jù)優(yōu)選實施例的方法中。將在隨后的節(jié)中說明關(guān)于該過程的三個步驟的細節(jié)。
“簡單算子分裂”的可選變型為 (i)以使用時間tn處的“舊”速度梯度
的關(guān)于Ω(n)的旋轉(zhuǎn)步驟開始，然后 (ii)進行“中間”步驟以將(i)的結(jié)果從Ω(n)擴展到Ω(n+1)，和 (iii)以使用時間tn+1處的流速U的對流步驟結(jié)束。
在這個變型中，以相反的順序執(zhí)行對流和旋轉(zhuǎn)步驟。盡管理論時間離散化誤差具有與上述第一種變型的情況相同的級，但是由于使用了不同瞬時時間的流速和速度梯度，因此這種可選變型趨向于一致性較差?？梢酝ㄟ^另一變型來避免該問題 (i)以如上所述的已知值ym(n)從Ω(n)到Ω(n+1)的擴展步驟開始，然后 (ii)執(zhí)行旋轉(zhuǎn)步驟，隨后 (iii)進行對流步驟 該變型對于步驟(ii)和(iii)使用Ω(n+1)中、時間tn+1處由流求解器提供的U和
的“新”值。
6.2“對稱算子分裂” 理論上來說，可以通過“對稱算子分裂”方法來獲得關(guān)于時間離散化誤差的更高階的精度。這種方法的基本想法是將具有完整步長的部分方程之一的一個中間步驟夾在具有半步長的其它方程的步驟之間。該方法的一種可能的變型導(dǎo)致以下的四步驟過程 1.擴展步驟從域Ω(n)的單元中的開始，計算在新域Ω(n+1)的所有單元中的初始值ym(n)。
2.預(yù)旋轉(zhuǎn)步驟從由擴展步驟提供的Ω(n+1)的已知值的擴展ym(n)開始，使用由流求解器提供的域Ω(n+1)中的速度梯度
通過具有半步長Δtn+1/2的ODE系統(tǒng)的數(shù)值解來計算預(yù)旋轉(zhuǎn)值ym(pre)。
3.對流步驟從由預(yù)旋轉(zhuǎn)步驟提供的作為初始值的ym(pre)開始，使用由流求解器計算的域Ω(n+1)中的流速U(rm，tn+1)，通過具有全步長Δtn+1的對流方程的數(shù)值解來計算中間值
4.后旋轉(zhuǎn)步驟從由前面的對流步驟提供的初始值
開始，使用由流求解器提供的域Ω(n+1)中的速度梯度
通過具有半步長Δtn+1/2的ODE系統(tǒng)的數(shù)值解來計算最終值 如上所述的“對稱分裂”變型產(chǎn)生完整方程的近似解，其具有步長的二階時間離散化誤差(即O(Δtn+12))。
另一種“對稱分裂”變型在理論上可能將旋轉(zhuǎn)步驟夾在具有半步長的兩個對流步驟之間。
(i)以將已知值ym(n)擴展到新域Ω(n+1)開始， (ii)使用流速U(rm，tn+1)在新域中進行半步長的“預(yù)對流”步驟，然后 (iii)使用速度梯度
在Ω(n+1)上執(zhí)行完整旋轉(zhuǎn)步驟，以及 (iv)以使用流速U(rm，tn+1)的另一個半步長“后對流”步驟結(jié)束，步驟(iv)最終導(dǎo)致也具有二階精度的的近似。如在“簡單分裂”的情況一樣，還存在第三種變型，其交換初始擴展步驟和預(yù)旋轉(zhuǎn)步驟，以后者開始并且從而使用Ω(n)上的速度梯度
的“舊”值。
由于在模具填充中遇到的典型流情況經(jīng)常包括導(dǎo)致纖維的顯著旋轉(zhuǎn)的剪切流，因此與純粹的對流輸運的影響相比，“旋轉(zhuǎn)算子”部分總是扮演重要(并且在大部分情形中為支配的)角色。因此通過針對第一對稱分裂變型詳細描述的預(yù)旋轉(zhuǎn)和后旋轉(zhuǎn)步驟的對流步驟的對稱夾入組成了“對稱分裂”方法的優(yōu)選變型[19]。然而，在大部分實際情況下，由流體流計算提供的時間步長的大小足夠小，并且6.1節(jié)所提出的“簡單分裂”變型在這些情況下以足夠的精度工作。
6.3可選分裂變型 基于右側(cè)函數(shù)F(...)的代數(shù)結(jié)構(gòu)，多種可選分裂方法是可能的。如果該函數(shù)由項的和F(...)＝∑kFk(...)組成，則可以將右側(cè)分裂成由部分函數(shù)Fk(...)確定的子方程的另外的集合。有可能某些部分函數(shù)包括線性算子，即在這種情況下，可以考慮部分對流方程(即)的右側(cè)的這些線性項，其由此保持線性，而將真正的非線性函數(shù)Fj(...)中的一些(或者全部)留在單個ODE或者如上所說明的單獨ODE的列表內(nèi)進行處理。
查看方程(16)的右側(cè)，可以看到三項的和，第一和第三項為線性的，而包括閉合近似的中間項可能具有關(guān)于FO矩陣的非線性依賴性。(基于閉合的選擇，可以將該中間項分裂成另外的子項和。)根據(jù)以上的評論，(16)的右側(cè)函數(shù)的項結(jié)構(gòu)開啟了構(gòu)建各種其它分裂方案的可能性。與基于分裂成方程(18a/b)的變型相比，這些其它替選方案不對應(yīng)于可相互分離的物理效應(yīng)(如“對流輸運”對局部“旋轉(zhuǎn)動力學(xué)”的情況中)的明確的區(qū)別。因此，它們僅僅可以被當(dāng)作對物理FO現(xiàn)象的適當(dāng)模擬而言幾乎沒有用的“人工數(shù)學(xué)分裂”。這里提及這些“人工”分裂變型的可能性僅僅是為了完整性并且在這里將不作進一步討論。
6.4“對流步驟”的數(shù)值方法 本節(jié)說明用于對對流方程(18a)進行積分的數(shù)值方法，需要將其作為閉合FTE的“算子分裂”方法中的子問題進行求解?；卺槍π碌臅r間步長計算的速度U(n+1)，使用一階逆風(fēng)方案[32]來組合輸運方程系統(tǒng)矩陣，并且針對域Ω(n+1)上的FO矩陣的每個分量求解所得到的方程系統(tǒng)。
6.5“擴展步驟”的數(shù)值方法 使用數(shù)值方法來將“舊”域Ω(n)的計算單元內(nèi)的FO矩陣的已知值擴展到“新”域Ω(n+1)的單元中。由經(jīng)修改的VOF方法執(zhí)行“新”域Ω(n+1)的計算，該VOF方法作為應(yīng)用在根據(jù)針對流方程(即質(zhì)量、動量和能量守恒方程)的數(shù)值解的優(yōu)選實施例的方法中的算法的一部分。對于兩個域共有的那些單元，在擴展步驟期間FO矩陣的值保持不變。一般而言，域Ω(n+1)包含若干“新填充的單元”，其需要初始分配FO矩陣的元素的合理值。根據(jù)優(yōu)選實施例，根據(jù)之前計算的質(zhì)量流，通過對應(yīng)于計算單元內(nèi)來自/去往其鄰近單元的質(zhì)量輸運的加權(quán)平均來完成該分配。
由于假定粒子濃度是均勻的，并且在宏觀水平上通過體積平均過程得到FO矩陣，因此假定單元將從其鄰近單元得到FO貢獻，該貢獻與輸運到其中的流體質(zhì)量的凈值成比例。根據(jù)之前總結(jié)的FTE的相空間(即所有可能的FO矩陣的集合)的數(shù)學(xué)特性，后者是被限制到(五維)超平面的實對稱3×3矩陣的六維向量空間的有界凸子集。由于根據(jù)質(zhì)量輸運加權(quán)的平均組成到FO矩陣的凸組合中，并且因此總是導(dǎo)致有效的FO矩陣，初始化過程與FT模型的相空間的拓撲結(jié)構(gòu)相容(理論上要求的)，這是該過程的重要性質(zhì)。
6.6“噴泉流”效應(yīng)的處理 術(shù)語“噴泉流”表征在一大類粘性、非牛頓流體的情況下的前沿處的自由表面附近的整體宏觀流型。在“噴泉流”中，流前沿附近的上行粒子從核心區(qū)域轉(zhuǎn)移到壁邊界。由于流體的材料性質(zhì)，這種效果實際上“自動地”發(fā)生并且不需要任何附加的建模，即至少在基于具有適當(dāng)?shù)姆桥ｎD本構(gòu)定律的Navier-Stokes方程的3D流計算中，“噴泉流”是突生現(xiàn)象。然而，如果流計算基于簡化模型(例如，如在[1]中的Crochet的文章的11.3節(jié)中所討論的Hele-Shaw類型)，則情況并非如此。在使用根據(jù)優(yōu)選實施例的方法的通道流(channel flow)模擬中，可以在粒子跡中清楚地辨認(rèn)“噴泉流”模式，這顯示由流求解器說明了這種現(xiàn)象。
針對新填充的單元的纖維取向的初始化過程，需要采取某種特殊的操作以確保垂直于自由表面的FO分量是零(與關(guān)于纖維取向的“噴泉流”的所觀察效果的一致性所要求的)。垂直FO分量的“被要求的”變?yōu)榱悴皇?有意地)由上述的“擴展步驟”過程所強制的，其原因如下從數(shù)學(xué)觀點來看，F(xiàn)TE模型是由雙曲輸運方程的右側(cè)耦合的雙曲輸送方程的系統(tǒng)(即具有非線性反應(yīng)項的對流反應(yīng)類型的系統(tǒng)，見上)。因此，不適合于(即數(shù)學(xué)上講不可能)將在自由表面垂直方向上的FO矩陣的分量歸零規(guī)定為針對自由表面單元處的FO矩陣的邊界條件。該FO分量的歸零應(yīng)該實際上是從所計算的流的“噴泉流”行為自動出現(xiàn)的。
然而，由于由時間和/或空間離散化誤差以及用于計算在模具填充過程中的自由表面的運動所引入的數(shù)值不精確性，有可能需要增加一些“校正”以確保針對新填充的單元的初始化過程與“噴泉流”現(xiàn)象所要求的自由表面處的FO行為一致。在所有新填充的單元上，檢查纖維取向張量的跡與理論值1相差不超過1％。如果差別太大，則通過正交投影到其特征空間中的特征向量的可容許的三角來校正纖維取向張量。
6.7模擬開始時的FO矩陣的初始化 (閉合)FTE的雙曲結(jié)構(gòu)要求在模具填充模擬的每個時間步、在澆口中(和附近)的計算單元中的FO矩陣的初始化作為邊界條件。因此需要做出這種初始FO狀態(tài)的合適選擇。在根據(jù)優(yōu)選實施例的方法中，出于以下所說明的目的選擇各向同性FO狀態(tài)(對應(yīng)于隨機FO分布)。
在高粘性剪切流中，F(xiàn)O狀態(tài)很大程度上受高剪切速率影響，因為它被快速地驅(qū)動到其局部準(zhǔn)平衡的“最終”狀態(tài)。因此在進?？?這是實際熔體進入該部分的地方)處所觀察的FO狀態(tài)由其貫穿澆注系統(tǒng)的流歷史完全確定，并且很大程度上不依賴于在澆口處假定的其初始狀態(tài)。另一方面，F(xiàn)TE的解析結(jié)構(gòu)及其運動學(xué)(即在相空間中可能的狀態(tài))的分析顯示，鑒于在澆口附近的澆注系統(tǒng)中的流場的隨后的適應(yīng)，澆口處的隨機取向的假設(shè)是最佳選擇，因為隨機取向狀態(tài)產(chǎn)生速度梯度的所有分量的全耦合。
上述算子分裂過程是參考圖11描述的簡化形式 步驟1設(shè)置新填充單元的初始條件和邊界條件， 步驟2根據(jù)流場移動纖維取向，和 步驟3根據(jù)速度梯度計算纖維取向。
7.“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)的數(shù)值方法 一種數(shù)值方法用來對ODE的耦合系統(tǒng)(18b)進行積分，需要將該耦合系統(tǒng)作為閉合FTE的“算子分裂”方法的一部分進行求解。已經(jīng)針對根據(jù)優(yōu)選實施例的方法的FO模塊的實施特別地開發(fā)了這一方法。它包含與短纖維增強熱塑性塑料材料的3D注塑成型模擬的典型的工業(yè)應(yīng)用情況下相關(guān)的各個方面，例如 ·使用FO矩陣的6個獨立元素的完整集的FTE的公式化(相對于使用跡條件以消除對角元素之一并將FTE的因變量的數(shù)量減少一個的“標(biāo)準(zhǔn)過程”)； ·在FTE的右側(cè)加入“懲罰”控制項，其將跡條件定義的超平面轉(zhuǎn)換為FTE的穩(wěn)定積分流形； ·利用與速度梯度分量相關(guān)的FTE的右側(cè)函數(shù)的特定縮放行為以構(gòu)造積分方法，該方法根據(jù)局部剪切速率的大小(速度梯度的最大分量)選擇時間積分方案； ·包含被限定在
區(qū)間內(nèi)的標(biāo)量取向因子的“穩(wěn)定混合閉合”的實施； ·使用最小數(shù)量的運算的用于具有混合閉合的FTE的特殊r.h.s.函數(shù)的估算的高效過程的實施。
下面詳細討論這些方面。
7.1混合閉合的穩(wěn)定化 由方程(6)和方程(6a/b)定義的混合閉合近似尤其在以下情況下遇到穩(wěn)定性問題速度梯度具有復(fù)雜的結(jié)構(gòu)，即不是簡單的剪切/拉伸類型而是反映復(fù)雜的3D流型，并且時間步長(如流求解器所確定的)的大小變得相當(dāng)大。標(biāo)量取向因子
向負值的偏離是這些不穩(wěn)定性的首要來源一旦該因子變?yōu)樨摚瑪?shù)值解很快變?yōu)椴环€(wěn)定，并且指數(shù)地偏離到遠離FO矩陣可容許的相空間集合的值。
因此對
的值的控制提高FO計算過程的數(shù)值穩(wěn)定性。
由于FO矩陣的行列式總被限制在
的區(qū)間內(nèi)，由方程(6c)定義的取向因子的理論上可容許的值也被限制。因此，以如下方式修改混合閉合的標(biāo)準(zhǔn)定義(6) (19) 由方程(6a/b)給出的線性閉合項和二次閉合項的定義保持不變，而由方程(6c)定義的取向因子
被上面的方程(19)中定義的受限取向因子
取代。
的該定義意味著如果
項被評估為區(qū)間
內(nèi)的值，則否則取0或1作為最小值或最大值。
定義(19)被稱為穩(wěn)定化混合閉合近似。數(shù)值試驗已經(jīng)顯示將
的值限制在區(qū)間
內(nèi)避免在所考慮的測試?yán)又挟a(chǎn)生嚴(yán)重的不穩(wěn)定。
已經(jīng)在各種3D注塑成型模擬中成功地測試了穩(wěn)定混合閉合。
7.2經(jīng)由跡條件減少變量的數(shù)量 作為方程(16)所示的閉合FTE的因變量的二階FO張量
為對稱3×3矩陣。由于其元素先驗地組成一組六個相互獨立的變量，因此FTE是六個PDE的耦合系統(tǒng)。雖然FTE的代數(shù)結(jié)構(gòu)形式上允許任意對稱3×3矩陣，但是該矩陣需要滿足表達為對其不變量的限制的若干附加條件以使其滿足作為FO矩陣的條件。
這些不變量條件中最簡單的跡條件開啟了明顯的可能性以消除FO矩陣的對角元素akk(2)之一，并且因此將變量的數(shù)量減少一個。大多數(shù)(如果不是全部)已發(fā)表的研究理論流情況(例如簡單剪切流或拉伸流)下的纖維懸浮液流的研究論文通過對于所要消除的特定對角元素作出固定的選擇來使用此方法。
然后通過將例如代入到方程(16)的右側(cè)來實現(xiàn)消除，該方程以這種方式變得不依賴于a33(2)。由于對于利用方程(14)或(15a)的所有合理的閉合近似FTE形式上與跡條件相容，因此跡條件和對角元素a33(2)的DPE一致地滿足(即作為形式上的代數(shù)恒等式)，并且只考慮FO矩陣的剩下的五個分量的集合以及進一步的它們的PDE的耦合系統(tǒng)就足夠了。此過程是純形式上的。對要以這種方式消除的對角元素的選擇是否是好的取決于流的局部速度場U及其梯度
的空間特性。對于僅顯示流速度的平面變化(例如在x1-x2平面內(nèi))的簡單類型的流，上面所說明的對a33(2)的消除是合理的選擇，因為在這種情況下纖維的主導(dǎo)旋轉(zhuǎn)動力被限制在流平面內(nèi)，因為在這種情況下速度梯度的特殊形式只導(dǎo)致對正交x3方向的“弱耦合”。根據(jù)這一理由似乎選擇消去a11(2)或a22(2)很可能會導(dǎo)致數(shù)值問題(例如關(guān)系到解的穩(wěn)定性)并且因此導(dǎo)致較不滿意的結(jié)果。然而，在復(fù)雜模具幾何形狀中的3D注塑成型模擬中情況則非常不同，因為流速U及其梯度
在空間和時間中以復(fù)雜的方式變化，因此不能假設(shè)其具有任何特別的簡單形式，所以固定的選擇(如上面的選擇)很可能不是最優(yōu)的選擇。因此通過跡條件消去FO矩陣的一個對角元素被認(rèn)為對于工業(yè)懸浮液流的模擬而言是不適合的。
在Shampine的一系列論文[34-36]中提出了反對通過單位跡條件消去一個對角元素的與數(shù)值穩(wěn)定性相關(guān)的另一個理由，在該系列論文中研究了具有特殊代數(shù)約束(“守恒定律”)的ODE系統(tǒng)。在[34-36]中考慮了ODE的耦合系統(tǒng)其中通過質(zhì)量守恒的條件將“物類濃度”的n維向量c＝(c1，...，cn)T限制在n-1維的超平面。而這個條件形式上產(chǎn)生了消去一個濃度的可能性，例如將關(guān)系代入右側(cè)函數(shù)F(...)，并由此同時精確地滿足守恒定律。Shampine指出在數(shù)值積分過程中使用這個“手法”僅僅導(dǎo)致濃度c1，...，cn-1的數(shù)值誤差的積累，其中濃度c1，...，cn-1通過(減少的)ODE系統(tǒng)的數(shù)值積分計算得出，其中經(jīng)由守恒定律代數(shù)地獲得剩下的濃度cn。雖然守恒定律由于構(gòu)建而總是被精確地滿足，但是無法保證以這種方式計算的數(shù)值解是準(zhǔn)確的。實際上，已知尤其對于剛性O(shè)DE系統(tǒng)來說，如果這種“直接消去方法”使用的不夠小心的話，數(shù)值解可能任意遠地偏離真實解。根據(jù)[34-36]，這些論據(jù)在包括權(quán)重因子的向量w＝(w1，...，wn)T的一般線性守恒定律的情況下同樣保持有效并且導(dǎo)致相同的結(jié)論，即便在考慮更為一般的具有g(shù)(c)＝const.形式的非線性限制的情況下也是如此。
由于通過“單位跡”限制(其為線性“守恒定律”)增強的Folgar-Tucker方程是[34-36]中出現(xiàn)的數(shù)學(xué)方程式的具體例子，Shampine的數(shù)學(xué)分析顯示，當(dāng)在FO矩陣的對角元素中引入不可控的偏差時，通過單位跡條件消去一個對角元素很可能導(dǎo)致數(shù)值不穩(wěn)定。
由于上面所給的理由，與通過單位跡條件消去一個對角元素的“標(biāo)準(zhǔn)方法”不同，根據(jù)優(yōu)選實施例的方法使用FO矩陣的全部六個元素以及FTE系統(tǒng)的方程用于FO計算。
7.3利用懲罰項的動態(tài)跡穩(wěn)定 如果閉合FTE系統(tǒng)及其“旋轉(zhuǎn)步驟”部分被認(rèn)為是由等式(16)和等式(18b)給出，并且不使用跡條件從獨立變量集合中消去FO矩陣的一個對角單元，則需要某些附加方法來保持?jǐn)?shù)值計算的FO矩陣的跡盡可能的接近其期望值1。
在論文[34-36]中建議避免任何類型的消去方法(如前面的7.2節(jié)所討論的)并且通過將通常的ODE積分方法與某些考慮(或精確或近似地)滿足所需要的守恒定律或代數(shù)限制的投影過程相結(jié)合來處理完整系統(tǒng)。對于每個時間步驟，首先在不考慮代數(shù)限制的情況下計算數(shù)值解，并且隨后通過向限制方程(例如線性限制情況下的超平面)所定義的代數(shù)流形的最接近點進行投影映射來進行校正。
作為此方法的替選的原則試圖通過向模型方程的右側(cè)加入適當(dāng)?shù)目刂祈椧宰鳛橄拗疲瑪?shù)值計算的MO矩陣的跡能夠被保持為盡可能接近其期望值1。控制項的存在在沒有任何附加措施的情況下產(chǎn)生限制或守恒定律的近似滿足。在許多情況下可以顯示出，“硬”控制有效地導(dǎo)致一種投影，而“軟”控制允許與所要求的限制間有較小的偏離。
在根據(jù)優(yōu)選實施例的方法中，控制方法尤其用于將數(shù)值計算的FO矩陣的跡近似地保持在其所需要的單位值。穩(wěn)定化是基于通過分別向等式(16)或等式(18b)的右側(cè)方程加入適合的控制項來將跡條件所定義的5D超平面變形到閉合FTE系統(tǒng)的穩(wěn)定積分流形的思想。從不同的角度來看，這等價于成為跡的微分方程(12)的穩(wěn)定固定點的情況。
已經(jīng)確定任何適當(dāng)?shù)目刂祈棏?yīng)滿足的一般要求如下 ·如果跡條件已經(jīng)被滿足，則控制項必須一致地變?yōu)榱恪?br>
·如果數(shù)值解不能滿足跡條件，則控制項必須迫使數(shù)值解回到超平面。
除了這兩個主要要求之外，合理的控制項應(yīng)該與和速度梯度的尺度改變(rescaling)相關(guān)的FTE方程的右側(cè)的縮放行為相容(關(guān)于這一點的更多細節(jié)參見關(guān)于流受控時間積分的7.4節(jié))，并且它應(yīng)該包含允許“軟”控制或“硬”控制的調(diào)節(jié)的調(diào)整參數(shù)，該“硬”控制實際上與“無限硬”控制的限制中的投影相似。
在混合閉合近似(標(biāo)準(zhǔn)和穩(wěn)定化變型)的情況下，對適合的控制項的特定選擇將適合于特定的閉合近似的選擇。從FTE開始，通過以下恒等式獲得跡的DE 在(標(biāo)準(zhǔn))混合閉合的特殊情況下，可以將右側(cè)函數(shù)的跡估算為
其產(chǎn)生跡的DE的簡單形式(15a)，其具有由(15b)給出的可變前因子
在穩(wěn)定化混合閉合的情況下，標(biāo)量取向因子
被等式(19)所定義的其受限的(穩(wěn)定化的)變型所代替，因此對于跡獲得了形式上完全相同的DE，其中前因子
包含
而非
在(穩(wěn)定化)混合閉合的特殊情況下對要加入右側(cè)函數(shù)
的懲罰項的具體選擇為
參數(shù)α需要為正，但不必是隨時間恒定的。3×3矩陣
需要為具有單位跡的對稱矩陣，其它方面是任意的。矩陣
定義控制項(20)的“方向”。它可以是隨時間恒定或可變的。后者包含將
定義為
的函數(shù)的可能。對矩陣
的不同選擇對應(yīng)于控制項的不同變型。根據(jù)優(yōu)選實施例的方法中使用的變型使用(即bij＝1/3δij)，其對應(yīng)于與超平面正交的方向上的控制項。一個可選的變型由以下選擇給出(見[18])。在加入控制項(20)之后，必須也考慮具有穩(wěn)定化混合閉合的閉合FTE的一般變型 相似地，控制項進入ODE(18b)的右側(cè)，其由此假設(shè)等式(21)的形式，其中實質(zhì)導(dǎo)數(shù)D/Dt由偏導(dǎo)數(shù)
代替。由于控制項(20)的跡估算為
因此對控制項的具體形式進行修改以根據(jù)下式的修改(15a/b) 當(dāng)α＞0時，值對應(yīng)于DE(22)的穩(wěn)定固定點，并且因此相應(yīng)的超平面變?yōu)榫哂谢旌祥]合的FTE的經(jīng)修改的(“控制的”)形式(21)的穩(wěn)定積分流形。作為結(jié)果，(21)的任意數(shù)值解的跡被控制項(20)保持在所要求的值1附近，這是因為該項迫使所有具有的解在
方向上趨向較低的跡值，并且相應(yīng)地迫使具有的解趨向較大的跡值(也在
方向)。
控制項(20)的“強度”是可以通過調(diào)整參數(shù)α而適當(dāng)?shù)恼{(diào)節(jié)的參數(shù)α的較小值導(dǎo)致相對“軟”的令跡值向1的校正，較大的α值導(dǎo)致較強的推(接近于瞬間投影的效果)向超平面。實際上，當(dāng)被顯式ODE積分方法使用時，較大的α值可能導(dǎo)致問題，因為在這種情況下較大的跡修正以交替的方式將數(shù)值解推向超平面的兩側(cè)，并因此導(dǎo)致偽振蕩。已經(jīng)顯示[18]在相對較寬的中間值范圍內(nèi)選擇α值在控制項不引入不希望的缺陷的情況下產(chǎn)生適當(dāng)?shù)臄?shù)值結(jié)果。
對于導(dǎo)致閉合FTE的相應(yīng)的右側(cè)函數(shù)
的閉合近似的任意具體選擇，控制項的一般定義 一致地導(dǎo)致跡DE(22)，并且因此導(dǎo)致對于(穩(wěn)定化)混合閉合的特殊情況的如上所述的超平面的穩(wěn)定化。如果我們允許對稱、單位跡3×3矩陣
為
的任意函數(shù)，則方程(20a)構(gòu)成控制項的最一般的形式。
7.4流受控時間積分 使用“算子分裂”方法，“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)(18b)的數(shù)值解需要在具有由流求解器提供的全局步長Δtn+1＝tn+1-tn的時間區(qū)間[tn，tn+1]中、在填充區(qū)域Ω(n+1)的所有計算單元(由空間向量rm標(biāo)記)內(nèi)計算。在每個單元中，局部速度梯度
的分量作為一組外部系數(shù)進入(18b)的右側(cè)。如果這些較大，則可以預(yù)期FO矩陣元的值在該時間區(qū)間內(nèi)變化顯著，而很小的速度梯度會導(dǎo)致FO矩陣元與由先前計算步驟提供的初始值相比可以忽略的變化。
在注塑成型模擬的所有實際情況中，在填充階段，速度梯度在模具型腔內(nèi)的填充區(qū)域顯著變化。已經(jīng)觀察到，剪切速率的值在接近腔壁處通常較高，而在壁間的空隙的中間區(qū)域處具有相當(dāng)?shù)偷闹怠＿@導(dǎo)致受到壁附近流速方向上較強剪切流的纖維較快地取向，而中心區(qū)域的纖維需要明顯較長的時間來改變其取向狀態(tài)，這導(dǎo)致眾所周知的“三明治結(jié)構(gòu)”(由壁附近的高取向區(qū)域間的具有相對較低纖維取向程度的“核”中心區(qū)域組成)，通常在由短纖維增強熱塑料制成的零件中觀察到該結(jié)構(gòu)。
用于ODE系統(tǒng)積分的數(shù)值方案通過適當(dāng)(優(yōu)選地適應(yīng)的)選擇步長來考慮右側(cè)函數(shù)的“強度”。試圖以這種方式避免由于以下原因造成的不準(zhǔn)確在較大值區(qū)間上的粗略分步或者右側(cè)函數(shù)的劇烈變化；以及在該區(qū)域上劃分過小的步長，其中右側(cè)函數(shù)的(絕對)值相當(dāng)小。通過適當(dāng)?shù)倪x擇步長以及積分方案，可以將ODE積分的數(shù)值誤差控制在期望的界限。
在“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)的特別情況下，不管由速度梯度的變化造成的右側(cè)方程如何改變，時間區(qū)間[tn，tn+1]上具有全局步長Δtn+1和全局?jǐn)?shù)值誤差εFO的數(shù)值積分必須對于計算區(qū)域的所有單元一致。為了這一目的，已經(jīng)設(shè)計出了特別的數(shù)值方案，其通過以下選擇以對右側(cè)函數(shù)的最少數(shù)量的估算實現(xiàn)這一目標(biāo) ·“積分規(guī)則”(歐拉、中點或四階Runge Kutta)以及 ·(無量綱)局部步長，以及 ·對于所選規(guī)則的此步長的局部步數(shù)量 其適合于由局部速度梯度值確定的局部右側(cè)函數(shù)的“強度”。本方案與用于數(shù)值ODE積分的“適應(yīng)步長控制”的已知方法不同，因為它利用與速度梯度分量相關(guān)的方程(16)或(18b)的右側(cè)函數(shù)的縮放行為并且其特別地適合于FTE系統(tǒng)。該方案固有的對于積分規(guī)則和局部時間步的數(shù)量和大小的選擇過程是基于所使用的積分規(guī)則的理論離散化誤差的試探方法。下面具體說明該方案的兩個方面。
方程(17b)或(18b)的右側(cè)函數(shù)
的具體形式可以從方程(16)獲得 為了清楚地解析對速度梯度
的所有依賴關(guān)系，需要代入有效速度梯度和剪切速率張量的定義(1a/b)，即 以及定義旋轉(zhuǎn)擴散參數(shù)和有效剪切速率的恒等式和公式對于如下考慮，抑制函數(shù)
經(jīng)由閉合函數(shù)
對FO矩陣的間接依賴性以及經(jīng)由有效剪切速率對速度梯度的間接依賴性以及對于恒定模型參數(shù)CI和λ的依賴性，由此將右側(cè)函數(shù)的代數(shù)表達重寫(如同已經(jīng)對方程(17c)的右側(cè)所做的)為以下的簡化形式 在
的變元表中只留下對FO矩陣
和速度梯度
的依賴性，這樣做將是有益的。使用右側(cè)函數(shù)(23)的ODE系統(tǒng)(18b)的等效變形(對應(yīng)于(17c))可以寫成以下形式 用某個因子σ＞0與
的分量的相乘
導(dǎo)致相應(yīng)的相乘
由于(23)右側(cè)的特殊的代數(shù)結(jié)構(gòu)，導(dǎo)致代數(shù)恒等式 表示函數(shù)
與其第二變元相關(guān)的縮放行為。
如前所述，一般任務(wù)是ODE的耦合系統(tǒng)(18c)在時間區(qū)間[tn，tn+1]內(nèi)的數(shù)值積分，其具有對于填充單元區(qū)域內(nèi)的所有單元一致的全局?jǐn)?shù)值誤差εFO。這可以通過利用右側(cè)函數(shù)
的特殊性質(zhì)(24)通過以下方式對速度梯度、右側(cè)函數(shù)和ODE系統(tǒng)(18c)進行縮放來實現(xiàn) ·以最大絕對值確定速度梯度分量的值，即(應(yīng)注意由于
，參數(shù)γmax還依賴于空間向量rm，其標(biāo)記局部計算單元；以及時間坐標(biāo)tn或tn+1，速度梯度在該時間坐標(biāo)進入右側(cè)函數(shù)。) ·引入縮放速度梯度以及縮放時間坐標(biāo)τ＝γmaxt(注意這兩個均為無量綱量，并且此縮放僅局部地應(yīng)用于單元rm以及時間區(qū)間[tn，tn+1]內(nèi))，并且使用經(jīng)縮放的量和等式(24)將ODE(18c)轉(zhuǎn)換為縮放形式 原始ODE(18c)需要在“全局”時間區(qū)間[tn，tn+1]上以“物理”步長Δtn+1＝tn+1-tn(例如以[s]為單位測量)積分，這對于區(qū)域Ω(n+1)的所有計算單元是相同的，并且右側(cè)函數(shù)的“強度”根據(jù)速度梯度
的變化而變化，如針對區(qū)域中的不同單元的參數(shù)γmax的不同值所指示的。與此不同，經(jīng)轉(zhuǎn)換的ODE(25)需要在具有無量綱大小的“局部”縮放的時間間隔[τn，τn+1]上積分 但是由于縮放的速度梯度
進入(25)的右側(cè)，由于下列事實，右側(cè)函數(shù)具有近似一致的“單位強度” i.縮放的速度
梯度由于構(gòu)建而為無量綱量，其分量的絕對值等于或小于1。
ii.FO矩陣元以及由閉合函數(shù)
提供的四階張量元素的值也被限制在1。
iii.因此當(dāng)用
估算時，右側(cè)函數(shù)
的所有分量都為O(1)階，而當(dāng)用
估算時
的最大分量具有與γmax相當(dāng)?shù)慕^對值。
使用任意矩陣范數(shù)‖...‖測量
的“強度”，可以將這些陳述給以準(zhǔn)確的數(shù)學(xué)形式 由于根據(jù)(27)的變形的ODE(25)的縮放的右側(cè)函數(shù)
對于區(qū)域Ω(n+1)內(nèi)的所有單元具有一致的“單位強度”，通過如下選擇可以在局部變化大小Δτn+1(m)的間隔上實現(xiàn)具有一致的全局誤差εFO的數(shù)值積分 i.具有足夠小的離散化誤差的“積分規(guī)則”，以及 ii.適合的子步數(shù)量以覆蓋經(jīng)縮放的時間間隔。
根據(jù)優(yōu)選實施例的方法中使用的積分方案使用一組屬于一步方法類(見[38]的16.1節(jié)和[39]的7.2.1-7.2.3)的簡單積分規(guī)則，通過引用將其包含于此。根據(jù)優(yōu)選實施例的方法的FO模塊中使用的具體積分規(guī)則為 ·簡單正向歐拉規(guī)則，其為一階準(zhǔn)確度方法并且每步需要右側(cè)函數(shù)的1次估算， ·“中點”或二階Runge-Kutta(RK2)規(guī)則，其為二階準(zhǔn)確度方法并且每步需要2次估算右側(cè)函數(shù)，以及 ·四階Runge-Kutta(RK4)規(guī)則，其為四階準(zhǔn)確度方法并且每步需要4次估算右側(cè)函數(shù)。
這些積分規(guī)則的依賴于其階次p的如下特性被用于構(gòu)建本發(fā)明的一個實施例中所用的特殊方案 a)p階的方法每子步需要p次估算右側(cè)函數(shù)。
b)如果以N個大小為h＝Δτ/N的等距子步在總大小為Δτ的間隔上對ODE(系統(tǒng))進行數(shù)值積分，則在該間隔的最終子步處積累的總誤差εtot可以被估計為εtot～Δτ·hp。
c)如果要求該總誤差小于預(yù)先選擇的閾值εmax，則子步的大小由h＜(εmax/Δτ)1/p限定。相似地，子步的數(shù)量需要大于N＞Δτ·(Δτ/εmax)1/p。
d)在總間隔上具有小于εmax的總誤差的積分可以通過用p階方法進行N個子步來執(zhí)行，其中間隔大小限制為Δτ＜(εmax·Np)1/(p+1)。這意味著如果滿足則大小為h＝Δτ的單子步(即N＝1)即足夠。
考慮到這些，可以構(gòu)建“混合”積分方案，該方案通過以下選擇在(26)中定義的大小為Δτn+1(m)的縮放的時間間隔上以小于εmax＝εFO的全局誤差產(chǎn)生方程(25)的數(shù)值積分 (i)具有三種積分規(guī)則之一的單步 (ii)使用RK4規(guī)則的適合大小的多個子步 其中該選擇取決于相比于所要求的誤差閾值εFO(見下面內(nèi)容)的Δτn+1(m)的相對大小。由于此誤差限制對于區(qū)域Ω(n+1)的所有計算單元是相同的，因此該積分以不依賴于局部值Δτn+1(m)(或相似的)的一致的誤差εFO執(zhí)行。
所提出的方法特定地選擇積分規(guī)則，該積分規(guī)則以根據(jù)上面的a)到d)方面所顯示的關(guān)系所估計的右側(cè)函數(shù)估算的最少次數(shù)來達到所要求的誤差限定。雖然由流確定的時間步長Δtn+1與總的注塑成型時間(其為秒量級)相比通常相當(dāng)小(Δtn+1≈10-2...10-4s)，無量綱量Δτn+1(m)可以為O(1)階，由于γmax的值可以因為小間隔尺寸以及聚合體熔體的高粘度而變得大得多(10...100s-1量級)。由于FO矩陣的更新值也為O(1)階，所要求的誤差界限的合理選擇在范圍內(nèi)改變，因此可以假設(shè)典型應(yīng)用中的子步長h≤0.1?？梢燥@示對于步長h＜1/2，RK4規(guī)則的步長h的單步比使用“中點”(RK2)規(guī)則的半步長h/2的兩步或四分之一步長h/4的四個歐拉步要更加準(zhǔn)確，雖然對于所有變型來說數(shù)值工作量(四次估算右側(cè)函數(shù))是相同的。
對積分規(guī)則和(子)步數(shù)量的選擇基于相對于針對值p∈{1，2，4}計算的一些列誤差界限的Δτn+1(m)的大小，該系列誤差界限定義各種積分規(guī)則的精度次序。由于0＜εFO＜1(實際上通常滿足εFO＜＜1)，誤差界限的大小隨著次序參數(shù)p的增大而單調(diào)增長(總滿足0＜εp＜εp+1＜1)，所以與下面給出的算法的公式化相關(guān)的誤差界限總是根據(jù)0＜εFO＜ε1＜ε2＜ε4＜1排序。(對于例如εFO＝0.001的典型的選擇，得到數(shù)值ε1≈0.032，ε2≈0.1以及ε4≈0.25。)相似地，對于無量綱縮放時間間隔Δτn+1(m)的大小有最小閾值εmin，在其之下即便步長為的單歐拉步僅產(chǎn)生FO矩陣元對于先前值的小到可以忽略的更新。這通常發(fā)生在速度梯度非常小的區(qū)域的所有單元中，例如在兩個相鄰壁之間相對較大的間隙尺寸的情況下的中心核區(qū)域中。例如εmin＝10-6的最小閾值對于典型模具填充應(yīng)用來說是合理的選擇。
考慮以上討論的所有方面，在區(qū)域的每個單元內(nèi)執(zhí)行的、并且在圖12的流程圖中圖解的算法可以通過如下方式進行表達 1.首先估算并且計算縮放的速度梯度和縮放的間隔大小 2.在步驟12.1中檢查是否 a.如果是，則跳過ODE積分(即什么也不做)，保持該單元的FO矩陣的先前值作為新(更新的)值并退出該過程。
b.如果否，則轉(zhuǎn)至步驟12.2。
3.檢查是否 a.如果是，則通過步長為的單歐拉步驟使用
估算右側(cè)函數(shù)來更新該單元中的FO矩陣的先前值，并退出該過程。
b.如果否，則轉(zhuǎn)至步驟12.3。
4.檢查是否 a.如果是，則通過步長為的單“中點”(KR2)步驟使用
更新該單元內(nèi)的FO矩陣的先前值，并退出該過程。
b.如果否，則轉(zhuǎn)至步驟12.4。
5.檢查是否 a.如果是，則通過步長為的單“四階Runge-Kutta”(RK4)步驟使用
更新該單元內(nèi)的FO矩陣的先前值，并退出該過程。
b.如果否，則轉(zhuǎn)至步驟12.5。
6.在的一般情況下，通過執(zhí)行具有適合步長的若干RK4步驟來執(zhí)行FO矩陣的更新。
a.達到所需要的εFO精度的最少步驟數(shù)量為滿足不等式N＞Δτ·(Δτ/εFO)1/p的最小整數(shù)N。使用函數(shù)int(...)，返回浮點數(shù)值的整數(shù)部分，以及函數(shù)max(...)，返回兩個數(shù)值中較大的一個，可以通過以下公式計算整數(shù)N≥1 N＝max(int(Δτ·(Δτ/εFO)1/p)+1，1). b.然后利用計算需要的步長。
c.一旦計算出N和h，通過RK4規(guī)則的具有步長h的N步(在右側(cè)使用
)更新該單元中的FO矩陣的先前值并退出該過程。
在一個優(yōu)選實施例中使用上述算法。它為“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)(18b)在區(qū)域Ωn+1中的所有單元在時間區(qū)間[tn，tn+1]上定義“混合”積分規(guī)則，其具有小于εFO的一致全局誤差和對(18b)的右側(cè)函數(shù)的最少估算次數(shù)。通過在右側(cè)函數(shù)的估算中使用縮放的速度梯度以及基于縮放的間隔長度Δτn+1(m)的大小計算步長h來“就地”完成該積分規(guī)則對(18b)的縮放形式(25)的應(yīng)用。
在此算法的典型應(yīng)用(具有誤差界限參數(shù)εFO＝0.001以及εmin＝10-6)中，觀察到在區(qū)域Ωn+1中的大多數(shù)單元(即大約80％)中，通過單歐拉步驟執(zhí)行FO矩陣的“旋轉(zhuǎn)步驟”更新，較少數(shù)量的單元僅僅由于(即γmax的較小值)而被“跨過”，并且有合理數(shù)量的據(jù)推測位于接近模具型腔壁并因此具有高剪切速率值的單元由一個或若干(通常不多于20個)“四階Runga-Kutta”步驟更新。
7.5具有動態(tài)跡穩(wěn)定化的“旋轉(zhuǎn)步驟”O(jiān)DE積分 如果右側(cè)函數(shù)包含控制項(混合閉合情況下的特殊類型(20)或者如方程(20a)給出的一般類型)以實現(xiàn)動態(tài)跡穩(wěn)定化，則右側(cè)函數(shù)的縮放行為仍然與本節(jié)描述的“混合”積分算法兼容，只要將方程(20)或(20a)的控制參數(shù)定義為α＝α0·γmax，其具有典型大小的無量綱控制參數(shù)α0～O(1)。
通過查看一般形式(20a)的控制項來對此進行說明。使用包含任意閉合近似的右側(cè)函數(shù)
的縮放特性(24)，我們可以根據(jù)等式的順序重寫它的跡 并且，使用控制因子α＝α0·γmax，得到控制項(20a)的等效表達式 右側(cè)的表達式現(xiàn)在包含作為單獨的因子的γmax以及依賴于
的項，該項由縮放的速度梯度
代替
進行估算，但是其它方面形式上與(20a)一致。這可以由以下等式表示 該等式還示出，如果將控制參數(shù)選擇為α＝α0·γmax，則在其一般形式(20a)下的控制項具有與利用γmax對速度梯度的縮放相關(guān)的“未經(jīng)控制的”右側(cè)函數(shù)
完全相同的縮放行為。
在(穩(wěn)定化的)混合閉合的特殊情況下，可以通過將項(20)重寫為以下形式來以同樣的方式從前因子
提取γmax
其中通過與
相同的公式使用
代替
計算“再縮放”的前因子
這明顯地顯示(28)對于在混合閉合情況下假設(shè)的控制項的特殊形式也是有效的。
以上的考慮顯示，針對“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)的“混合”積分算法在不被改變的情況下也能夠應(yīng)用于有控制項的情況，只要α0用作控制參數(shù)并且用縮放的速度梯度估算依賴于右側(cè)函數(shù)的項。簡單地通過向方程(25)右側(cè)加入縮放的控制項來得到對應(yīng)于方程(21)的縮放形式的ODE，即需要考慮用ODE 代替(25)作為積分算法的分離時間步驟的基礎(chǔ)。
在根據(jù)優(yōu)選實施例的方法的FO模塊中實施的時間積分方案通過將本節(jié)中定義的“混合”積分算法應(yīng)用到包括利用控制項(20)進行的動態(tài)跡穩(wěn)定化的ODE系統(tǒng)(29)來執(zhí)行“旋轉(zhuǎn)步驟”O(jiān)DE系統(tǒng)的數(shù)值積分。
7.6具有穩(wěn)定化混合閉合的FTE的高效估算 右側(cè)函數(shù)的估算是FO計算過程中消耗最多的部分，因此本方法以最高效的方式進行估算。
通過7.4節(jié)提出的“流受控時間積分”(FCTI)的算法方法已經(jīng)部分地解決效率方面。使用該方法，能夠以對FTE右側(cè)函數(shù)(可能包含控制項)的最少估算達到一致精度。
直接影響FO計算的計算消耗的最重要的因素是對閉合近似的選擇。穩(wěn)定化混合閉合以低計算消耗產(chǎn)生合理的精度，并且被本方法采用。
具有混合閉合的FTE的代數(shù)變形 提高效率的第一步是將混合閉合的代數(shù)定義(6)代入FTE的右側(cè)函數(shù)(23)，(在幾個代數(shù)變形之后)產(chǎn)生特定代數(shù)形式的右側(cè)函數(shù)
(30)的右側(cè)的前兩項的和在形式上與表達式
一致，但是在1(a)給出的有效速度梯度矩陣
處包含矩陣 右側(cè)剩下的三項都具有標(biāo)量前因子
與矩陣量
的乘積
的形式。該前因子由一組公式給出
以上給出的針對前因子
的公式與(15b)中給出的公式一致，因為并且其中 公式(32)在可壓縮或不可壓縮的流速場的情況下有效。這里強調(diào)，雖然理論上當(dāng)假設(shè)流為不可壓縮時應(yīng)滿足divU＝0，這不被明確地使用，而是將所有包含的項保持在公式(32)中。如果由于數(shù)值誤差(其不可避免地發(fā)生在流計算過程中)導(dǎo)致發(fā)現(xiàn)(32)將導(dǎo)致穩(wěn)定行為，而當(dāng)在前因子公式(32)中將與
成比例的項先驗地假設(shè)為零的情況下則會發(fā)生不穩(wěn)定。
形式(30)下的具有量(31)和(32)的右側(cè)函數(shù)的使用是通向高效估算的基礎(chǔ)步驟，因為右側(cè)函數(shù)(30)的代數(shù)結(jié)構(gòu)被以這樣的方式組織使得各種計算只需要進行一次(當(dāng)計算前因子時)，這節(jié)省了很多運算。
若干代數(shù)運算發(fā)生在頻繁出現(xiàn)的共同的子表達式(例如(31)和(32)中的項2Dr或(1-fs)/7)中，只對它們進行一次計算并將它們的值存儲在輔助變量中用于后面的使用，這顯著地減少了計算工作量。通過接下來的步驟，即根據(jù)
重新定義矩陣
并且使用矩陣
代替
進行右側(cè)方程求解，省去了另外的一些運算。利用代數(shù)恒等式
省略了項
的詳盡計算，并且將右側(cè)函數(shù)根據(jù)下式重新定義
利用方程(31a)計算
的附加運算為(i)1個乘法以得到
以及(ii)3個“加法”以從矩陣
的對角元減去
所節(jié)省的運算數(shù)量等于FO矩陣的所有分量與前因子
的相乘以及
的減法運算。因此，如果使用方程(30a)和(31a)代替(30)和(31)用于估算右側(cè)函數(shù)，則運算的總數(shù)量更少。
在穩(wěn)定化混合閉合的情況下，根據(jù)(19)計算受限變量取向因子
因此，用
代替(32)中的fs在不影響計算消耗的情況下得到“穩(wěn)定化前因子”
和
控制項(見等式(20))
的加入導(dǎo)致對于右側(cè)估算的另外的運算需要。此控制項的一般結(jié)構(gòu)由具有前因子
的
給出，并且因此具有如上面討論的相同的形式
如果為投影矩陣選擇或則前因子
被并入前因子
或
之一。在第一種情況下，項
被加入到前因子
導(dǎo)致經(jīng)修改的右側(cè)函數(shù)
第二選擇導(dǎo)致需要加入到前因子
的附加項
右側(cè)估算的高效算法 利用上述步驟，顯示出一種方案，其根據(jù)以下計算序列以最小數(shù)量的代數(shù)運算估算(30b) 1.輸入
和參數(shù)λ、CI、α 輔助變量
λ+，λ-，fs，
φ1，ζ，ξ1，...，ξ6 2.由利用(1a)計算的有效速度梯度矩陣初始化
λ±＝1/2(λ±1) 3.計算對稱矩陣和它的跡 4.根據(jù)壓縮由計算ξ2＝2Dr。
5.計算標(biāo)量取向因子其括號形式以及項 6.計算壓縮 7.通過以下公式使用存儲在輔助變量ξk中的值根據(jù)(32)計算前因子
8.計算項
然后計算
9.使用輔助變量ξ5←2ξ3和
計算矩陣
應(yīng)該通過以下運算序列進行計算 10.通過計算初始化右側(cè)結(jié)果矩陣。然后相繼完成右側(cè)函數(shù)(30b)的估算
11.結(jié)果 使用“簡化符號”(CN)的向量形式FTE的重新表示。
右側(cè)函數(shù)的高效估算過程的構(gòu)建中的另一個要素基于以下事實FO矩陣
和右側(cè)函數(shù)
都是對稱矩陣并且只有6個獨立矩陣元。因此將它們作為一般3×3矩陣存儲是不實用的，并且將估算
所需要執(zhí)行的運算作為對矩陣的運算來實現(xiàn)是低效的。根據(jù)本發(fā)明，通過使用以下識別方案將對稱3×3矩陣作為實6分量向量進行處理向量索引μ∈{1，...，6}，以及對稱矩陣的索引對(ij)＝(ji)，其被稱為“簡化符號”(CN，參見例如[22])(在結(jié)構(gòu)力學(xué)中，這也被稱為“Voigt符號”) 這個方案的特別選擇當(dāng)然只是6?。?20個等效方案中的一個。上面的表格顯示通常的選擇，這也在優(yōu)選實施例的FO模塊中采用。(可替選的變型通常選擇非對角矩陣元對向量索引μ∈{4，5，6}的不同的分配。)使用CN方案產(chǎn)生對稱3×3矩陣的元素到相應(yīng)的6矢量(即
的向量)的分量的雙射映射
如下面所示
這樣，存儲右側(cè)函數(shù)
估算結(jié)果的FO矩陣、矩陣
和矩陣
被分配到相應(yīng)的向量。以同樣的方式，如果四階張量
關(guān)于第一和第二索引對對稱的話(即，如果滿足)，則該張量的分量被分配到6×6矩陣的元素
如果除此之外該張量還關(guān)于兩個索引對的互換對稱并且因此具有正交性對稱(見方程(11))，則滿足
即矩陣
本身對稱并具有個獨立元素。通過進一步的對稱性和張量分量間的代數(shù)關(guān)系，可以減少獨立元素的數(shù)量。在全對稱四階張量的情況下，存在6個附加恒等式
其將
的獨立元素的數(shù)量減少到15。規(guī)范化條件和跡條件產(chǎn)生一組代數(shù)關(guān)系，將獨立矩陣元的數(shù)量減小到更小的值(詳見[22])。
引入CN方法導(dǎo)致如下的“向量形式”的FTE
對于任意矩陣3×3矩陣
映射
定義實對稱3×3矩陣的六維向量空間中的線性映射，其可以被形式上寫成由實6×6矩陣
描述的
中的矩陣矢量乘積
其元素或者為零或者由
的元素的簡單代數(shù)項給出。
在這個意義上，CN方案產(chǎn)生標(biāo)識
其中
相似地，由以下公式以元素方式定義壓縮運算
如CN方案所給出的，使用索引分配此運算還可以寫成矩陣—向量乘積
其中
的矩陣元與
的矩陣元經(jīng)由如下形式相關(guān)
這揭示矩陣
(與
不同)不是對稱的！根據(jù)CN方案(即
)，矩陣
通過分配張量
的元素而產(chǎn)生，該張量本身通過閉合近似從FO矩陣計算得出，因此
項代表CN中的閉合近似。
所提出的算法的計算消耗 針對右側(cè)函數(shù)的高效估算而提出的過程的步驟2到步驟10中所執(zhí)行的計算操作的最高效的實施通過將CN方案應(yīng)用到此過程而產(chǎn)生。下表給出執(zhí)行使用CN方法的過程的各個步驟所需要的運算(
乘法和除法的數(shù)量，
加法和減法的數(shù)量)的數(shù)量和函數(shù)調(diào)用的概況
以下要點評述實施的各方面以及針對上述算法的各步驟在表中給出的運算計數(shù) ·由于矩陣
不是對稱的，因此將其作為3×3矩陣存儲。如果以
的計算初始化
則不需要
的額外存儲。
·用于
的計算，對稱矩陣
需要為矩陣形式(→步驟9)，但是在最終步驟10中，其還作為6向量出現(xiàn)。
·存儲和分配運算沒有被計數(shù)，因為它們對總的計算消耗的貢獻是可以忽略的。
·通過重新使用已經(jīng)存在的變量(當(dāng)不再需要其值時)來減少輔助變量的數(shù)量。這沒有在上面提出的算法中執(zhí)行，因為這會降低對算法進行的說明的清楚程度。
·對一對對稱矩陣的壓縮運算(如步驟4和步驟6中所做的)是使用CN通過7個
運算和5個
運算實現(xiàn)的。(在矩陣符號中，它由定義。) ·實對稱3×3矩陣的行列式(→步驟5)是使用CN通過11個
運算和5個
運算計算的。
·如果使用輔助變量將
作為對稱3×3矩陣進行存儲，則步驟9中執(zhí)行的矩陣
的組合需要6個
運算和12個
運算，外加2個乘法以計算變量ξ5和ξ6。
·矩陣運算需要使用CN的27個
運算和21個
運算。
·任何加(或減)多個單位矩陣的
類型的運算只需要β到對角元的三個加法。
表的最后一行示出獲得右側(cè)函數(shù)的單次估算需要96個
運算、84個
運算和3個函數(shù)調(diào)用的計算消耗以計算(雙精度)浮點數(shù)值的平方根以及一對這種數(shù)值的最小值和最大值。在當(dāng)前計算機硬件上，加法和乘法運算的計算消耗大致相同，對min或max函數(shù)的調(diào)用消耗大約1.5-2次運算，而平方根計算的消耗總共約為6-10次運算(即，平均8次運算，取決于所需要的精度)。
利用以上提出的算法實現(xiàn)的(30b)的單次估算總共需要190次運算的計算消耗。這是估算具有包括動態(tài)跡穩(wěn)定化的(穩(wěn)定化)混合閉合的FTE右側(cè)的最高效的方法。
計算效率評估 根據(jù)與估算由方程(35)給出的向量形式的右側(cè)函數(shù)的“標(biāo)準(zhǔn)”方法進行的比較產(chǎn)生對所提出過程的計算效率的評估。在以下情況下會需要采用這種方法，例如，如果編碼被規(guī)定為具有特定模塊結(jié)構(gòu)，使得(i)以“一般”過程執(zhí)行不依賴于閉合近似的所有運算，這種過程將6×6實矩陣
作為輸入，并且(ii)使用如上面所說明的CN方案和(35b)根據(jù)閉合函數(shù)
的具體選擇以單獨的過程計算
的矩陣元?？梢酝ㄟ^如下考慮來估算這種“標(biāo)準(zhǔn)”方法的計算消耗 ·該“標(biāo)準(zhǔn)”過程將從計算矩陣
6矢量
和標(biāo)量參數(shù)2Dr(見步驟2到步驟4)開始，加上平方根計算，這總共需要46次運算。
·如果事先解析地完成矩陣向量乘積并且只對所產(chǎn)生的針對六個分量的公式進行數(shù)值估算，則對應(yīng)于
的運算
需要48次運算。
·矩陣量乘積
以及將其從右側(cè)減去的計算需要66+6＝72次運算。
·2Dr·1的加法和6Dr·a的減法(使用6Dr＝3/2·2Dr)需要另外3+12+1＝16次運算。
為平方根計算計8次運算，到這時該方法用于估算(35)右側(cè)的計算消耗合計達190次運算。用于跡穩(wěn)定化的控制項的加入增加了更多的運算(見步驟8)，但是通過將項3Dr＝3/2·2Dr從
的對角線減去以代替從右側(cè)直接減去6Dr·a可以節(jié)省大約同樣多的運算。
這些“固定消耗”不依賴于計算矩陣
的元素所需要的閉合近似，該計算需要在單獨的過程中完成。
由于(6)和(6a/b/c)所定義的張量函數(shù)
具有正交對稱性，但由于二次項(6a)而不是完全對稱的，因此在混合閉合的情況下矩陣
有21個獨立元素。將
的定義轉(zhuǎn)換到CN產(chǎn)生計算矩陣元
時需要使用的如下一組表達式
對稱矩陣
和
都包含需要事先解析計算的全對稱四階張量表達式的分量。
當(dāng)矩陣元
為常量時，矩陣元
為向量分量aμ的線性表達式。矩陣元
和
的解析計算產(chǎn)生如下一對實對稱6×6矩陣
使用這些矩陣，根據(jù)如下過程確定組合矩陣
的計算消耗 ·由有效算法的步驟5完成受限標(biāo)量取向因子
和輔助變量ζ2的計算并且需要23次運算(為內(nèi)嵌的min(0，max(...，1))調(diào)用算作3次運算)。由一次附加除法計算輔助變量ζ1＝ζ2/5，因此計算
和ζ1/2所需要的運算量為24。
·由
初始化矩陣
這對于上三角的21個獨立矩陣元中的每個進行2次乘法，或總共42次運算。(利用對稱性分配剩下的矩陣元。) ·
的矩陣元由單次運算(即3ζ1的計算)以及將ζ1或3ζ1到相應(yīng)的非零矩陣元的若干分配來計算。相似地，整個更新運算
只需要從
的相應(yīng)條目減去ζ1或3ζ1。由于只需要在
的上三角元素執(zhí)行這些減法，因此更新運算的總運算計數(shù)降到10次(即9次減法加上1次計算3ζ1的運算)。
·矩陣
包含12個不同的表達式，其中只有9個需要由單次加法或乘法計算。(剩下的通過分配獲得。)由于12個表達式都不能被先驗地假設(shè)為零，
的計算需要12次乘法，因此對于該矩陣的獨立元素，計算
的總運算量為21次。隨后通過分配得到剩下的矩陣元。
·
的組合過程中的最終步驟由更新運算
組成，其需要21次運算以將矩陣
和
的上三角相加并進行分配以獲得剩下的下三角矩陣元。
根據(jù)(36)組合矩陣
的完整過程按順序執(zhí)行以上所描述的步驟，其需要24+42+10+21+21＝118次運算。(應(yīng)該注意到此過程也是以非常高效的方式建立的！)由于非對稱矩陣
是根據(jù)(35b)利用18個乘法從對稱矩陣
計算得出，因此用于上述計算矩陣
的過程的總運算量為118+18＝136次。
當(dāng)應(yīng)用到混合閉合式時，以上所述的“標(biāo)準(zhǔn)過程”總共需要約326次運算。這達到所提出的高效過程的計算消耗的大約1.72倍(即超出大約70％)。使用高效過程以“標(biāo)準(zhǔn)”方法的計算消耗的大約60％產(chǎn)生對于具有穩(wěn)定化混合閉合(包括動態(tài)跡穩(wěn)定化)的FTE的右側(cè)函數(shù)的單次估算的結(jié)果。
計算效率的這一提高主要是利用基于先前
的解析計算來巧妙地重組右側(cè)函數(shù)的具體代數(shù)結(jié)構(gòu)而達到的，其中通過經(jīng)由CN方案減少方程來最終達到優(yōu)化的效率水平。
7.7跡尺度改變、不變量監(jiān)測和相空間投影 任何FO矩陣需要具有非負特征值和單位跡(→2.4節(jié))。這意味著對于FTE的解的一組代數(shù)限制，該FTE由此變?yōu)镈AS(→2.6節(jié))。除了針對跡的等式針對特征值的非負性限制可以等效地以針對矩陣的第二和第三不變量和的一對不等式Ka≥0和Da≥0(→方程(9)、(10))的形式進行公式化。
通過加入FTE右側(cè)的適合的懲罰項經(jīng)由動態(tài)跡穩(wěn)定化(DTS)以“主動的”和計算上花費不多的方式將跡條件計算在內(nèi)，其中該懲罰項自動地將數(shù)值解的跡保持在所需要的值Sa＝1附近。使用此技術(shù)僅獲得近似地滿足跡條件(即Sa≈1)的FTE的數(shù)值解。從實際的觀點看，這幾乎不影響解的正確性，因為總有可能通過乘以因子1/Sa來改變矩陣元的尺度。將這種跡尺度改變運算數(shù)學(xué)地描述為
并具有若干有利特性 ·如果μk為矩陣
的特征值，改變尺度的矩陣
的特征值由μ′k＝μk/Sa給出，因此(i)如果Sa＞0，則它們不改變符號，并且(ii)如果Sa≈1(如DTS所保證的)，則它們的數(shù)值僅被輕微縮放。
·相應(yīng)的特征向量不被尺度改變運算影響。
因此跡尺度改變的運算不導(dǎo)致計算的數(shù)值解的質(zhì)的改變，而是在保留矩陣的所有本質(zhì)特性的情況下對解進行輕微修正。從這個觀點來看，F(xiàn)TE的具有非負特征值μk(或等效地非負的第二和第三不變量Ka和Da)且Sa≈1的任意數(shù)值解在實踐中可以被看做是FO矩陣。
FTE的不變量控制的一般方面 根據(jù)[12]的結(jié)果(→第3節(jié))，如果實對稱3×3矩陣的跡為正并且第二和第三不變量Ka和Da為非負，則保證了該矩陣特征值的非負性。
由于矩陣的跡是矩陣元素的簡單線性函數(shù)，并且實質(zhì)導(dǎo)數(shù)算子與跡運算對易，因此有可能以直接的方式得到跡的演化方程(→第4節(jié))。在“精確”情況以及混合或一般正交閉合近似的情況下(見方程(14)和(15a/b))，由于FTE右側(cè)的特殊代數(shù)結(jié)構(gòu)，此演化方程采取閉合形式。由于FTE的這些特殊代數(shù)性質(zhì)，有可能向FTE右側(cè)本身加入懲罰項，其以特殊(預(yù)期的)方式(→DTS)影響跡的演化方程。
然而，與跡不同，第二和第三不變量Ka和Da是矩陣元的非線性函數(shù)。因此不可能應(yīng)用相同(或相似)的過程來獲得針對這些不變量的閉合形式演化方程(通常不變量的演化方程的右側(cè)是所有矩陣元的函數(shù)，不只是其不變的組合，并且不獨立于編碼于特征向量中的矩陣方向特性)。相似地，很難(如果不是不可能的話)在FTE右側(cè)加入以預(yù)定方式影響非線性不變量演化的任何特定項。因此，通過FTE右側(cè)適當(dāng)?shù)膽土P項以將不變量保持在非負區(qū)域內(nèi)的第二和第三不變量的“主動”控制是不可用的。
一個可替選過程由“不受控”(或部分受控)積分過程與將數(shù)值解映射回其容許區(qū)域的投影運算的結(jié)合組成。這種“被動”過程導(dǎo)致有效的積分方案，如果數(shù)值積分的運算和隨后的投影被應(yīng)用在每個時間步，則該積分方案保持“不受控”方案的精度(具體討論參見[9]的IV.4章關(guān)于投影方法的部分)。
應(yīng)用到FTE，此結(jié)合的過程由以下方面組成 1.使用第5、6和7節(jié)(直到前面的7.6小節(jié))所描述的方法進行的積分步驟，然后 2.不變量監(jiān)測步驟，用來確定在先前的積分步驟中計算的矩陣是否仍為FO矩陣，以及 3.相空間投影步驟，如果監(jiān)測步驟的結(jié)果為否定(即違反了不變量條件)，則將矩陣映射回相空間集MFT。
由于相空間MFT為凸集(→第3節(jié)中的定理1)，對于每個實對稱3×3矩陣，MFT中存在“最小距離”(以適合的度量測量，見下面)上的唯一的矩陣。這定義了上述過程的第三步驟所需要的投影映射。下面討論相空間投影的更多細節(jié)。
高效不變量監(jiān)測和特征值計算 根據(jù)上述過程，需要在每個時間步和計算區(qū)域的每個(填充)單元中執(zhí)行不變量監(jiān)測。因此不變量(第二和第三不變量或特征值，其本身也是不變量)的計算需要以最高效的方式實現(xiàn)以避免任何不必要的計算開銷。使用符號(x，y，z)和(u，v，w)來表示第3節(jié)中引入的矩陣
的對角元和非對角元，由下列公式將不變量作為矩陣元的多項式表達式給出 其總共可以用11次加法和13次乘法估算。由于在FTE的數(shù)值積分過程中的每個時間步都要計算跡Sa，因此可以對于每個單元以每時間步另外的22次運算的消耗來完成通過檢查不等式Ka≥0和Da≥0進行的不變量監(jiān)測。
對不變量Ka和Da的檢查是用來確定實對稱3×3矩陣是否屬于FO矩陣的MFT集的最高效的方式，因為以條件μk≥0對特征值進行的檢查需要對后者進行詳盡的計算。
具體分析顯示可以通過計算特征多項式Pa(μ)的根來最高效地完成矩陣
的特征值的數(shù)值計算，其在計算作為Pa(μ)的系數(shù)的不變量所需要的24次運算之外消耗大約100次運算，因此其消耗高出大約5倍。這量化了使用特征值代替使用不變量用于監(jiān)測FO矩陣特性的計算開銷，并且突出了使用不變量(38)是高效監(jiān)測
的譜性質(zhì)所選擇的方法。矩陣的數(shù)值對角化是最昂貴的計算特征值的方法。因此不推薦這種過程，除非也關(guān)心特征向量(例如作為下面討論的相空間投影過程的一部分)。
相空間投影 如果不變量監(jiān)測指示在當(dāng)前時間步中FTE的數(shù)值解移出了相空間集MFT(即違反不等式Ka≥0和Da≥0中的一個或兩個)，則需要通過向相空間MFT的投影校正該解。
對于任意3×3實矩陣
通過定義“最接近”矩陣來產(chǎn)生該投影映射，其中以Frobenius規(guī)范測量兩個3×3實矩陣的距離其中Frobenius規(guī)范由3×3實矩陣的向量空間中的標(biāo)量積自身導(dǎo)致。這些量對對稱矩陣的六維子空間的約束是直接的?？梢燥@示出[18]對于任意3×3實矩陣
存在唯一矩陣使得
和MFT的元素
之間的距離
準(zhǔn)確地在的情況下最小化。
由矩陣給出該最小距離解，該矩陣由以下特性唯一地定義(其形式證明參見[18]) ·
的特征向量與矩陣
的特征向量相同。
·
的三個特征值(μ1*，μ2*，μ3*)是三角集(見第3節(jié)中的等式(7)和圖7)
的唯一點，該點具有
中距離由矩陣
的三個特征值給出的點(μ1，μ2，μ3)最小的歐幾里德距離。
在對稱3×3矩陣的六維實矢量空間
上定義的并且具有MFT中的值的投影映射
由如下公式給出簡潔的數(shù)學(xué)符號
并且實際計算投影映射的值的算法由以下步驟組成 1.給定實對稱3×3矩陣
作為輸入，計算其不變量Sa、Ka和Da。
2.檢查條件Sa＝1、Ka≥0和Da≥0。如果所有條件都滿足，則已經(jīng)成立，并且投影映射僅是等式
在這種情況下將輸入矩陣
分配給輸出并退出。否則(即如果Ka＜0或Da＜0)執(zhí)行下列步驟 3.計算該矩陣的特征值μk和相應(yīng)的特征向量Ek(即)。
此步驟可以形式地表示為 4.然后尋找對于
中的μ＝(μ1，μ2，μ3)具有最小歐幾里德距離的唯一點即 5.最后通過以下公式組合投影映射的值
注意到如果則(39c)產(chǎn)生μ*＝μ，因此在這種情況下(39d)符合(39a)。在這個意義上講，(39b/c/d)已經(jīng)表示了投影映射的一般定義。僅對于的情況需要執(zhí)行計算(39b/c/d)的結(jié)果的各個步驟，而在的情況下(39a)使該定義形式上完整。
雖然上面給出的計算
的算法的描述是數(shù)學(xué)上精確的，但是該算法的實際實施取跡Sa≈1的輸入矩陣(如具有DTS的FTE時間積分方案所提供的)，這將略去本過程的步驟1中對跡的附加估算。另外，實際實施將并入各個算法步驟的下列變化和/或說明，見圖13 ·在步驟1中，由(38)計算輸入矩陣的不變量Ka和Da。
·在步驟2中，只檢查條件Ka≥0和Da≥0，并且如果兩個條件都滿足，則使輸入矩陣保持不變并將其分配到輸出。
·通過數(shù)值矩陣對角化獲得步驟3的結(jié)果，例如通過迭代循環(huán)雅可比旋轉(zhuǎn)(見[38]，11.1章或[39]，6.5.2章)或單Givens/Housholder約簡步驟以及隨后的QR迭代(見[42]或[38]的11.2章)。
·實際上只在至少一個輸入特征值μk為負時執(zhí)行步驟4，所以求解點總位于三角集Δa的邊緣(包括角點)。
使用“跡尺度改變的”特征值(見方程(37)和隨后關(guān)于此問題的敘述)，將根據(jù)(39c)確定μ*所需要解決的最小化問題限制到由三角形Δa的三個角點(1，0，0)、(0，1，0)和(0，0，1)固定的平面Sa＝1。可以在
中等效地解決此“平面”最小化問題，例如在(μ1，μ2)平面中，通過尋找位于投影三角形邊緣上的唯一的對(μ1*，μ2*)
(即Δa到(μ1，μ2)平面的正交投影)，該唯一的對具有到對
的最小的歐幾里德距離。μ*的第三坐標(biāo)因此由給出。因為通過DTS，Sa≈1，解決“平面”最小化問題的算法產(chǎn)生對(39c)的精確解的很好的近似，而這比在
中直接根據(jù)(39c)計算μ*的算法簡單的多并且能夠更高效地實施。
·在步驟5中，計算矩陣
的各個元素的公式由給出，其以旋轉(zhuǎn)矩陣的矩陣元Rij的方式給出，該矩陣元已經(jīng)在步驟3中事先算出。
8.實施 該過程的結(jié)果為描述物件的給定區(qū)域中纖維取向可能性的分布函數(shù)，該結(jié)果以圖形或數(shù)字的形式顯示在計算機工作站(未示出)的顯示器上。這可以是執(zhí)行計算的計算機或工作站的顯示，或者可以在與執(zhí)行模擬的計算機網(wǎng)絡(luò)連接的計算機的顯示器上。
模具或產(chǎn)品的設(shè)計者將使用該模擬的結(jié)果來改進由注塑成型過程得到的物件的質(zhì)量。也可以由工程師利用該模擬的結(jié)果來降低用于制造所考慮的物件的成本。關(guān)于纖維信息的可靠信息的優(yōu)點使得工程師能夠在確定該物件的強度、剛性和穩(wěn)定性特性之前得到更好的理解和信息，因為纖維(纖維通常比聚合物材料具有更高的強度特性)的取向?qū)τ谖锛膹姸?、剛性和穩(wěn)定性特性具有決定性的影響。另外，纖維取向?qū)τ趯ζ渲袘腋∮欣w維的聚合物團進行注塑成型時可能發(fā)生的翹曲效應(yīng)有影響。通過得知纖維取向，造成影響的翹曲和其它應(yīng)力能夠被更好地預(yù)測或者設(shè)計改變以避免。
該模擬的結(jié)果還可以被傳輸?shù)狡渌鼞?yīng)用軟件，例如CAD軟件。這樣，模擬的結(jié)果可以用于設(shè)計軟件和模擬軟件之間的交互式過程。
本發(fā)明具有很多優(yōu)點。不同的實施例或?qū)嵤┛梢援a(chǎn)生一個或更多如下優(yōu)點。應(yīng)該注意到，這不是完全的列舉，可以有此處未描述的其它優(yōu)點。本發(fā)明的一個優(yōu)點是可以以顯著減小的計算工作量確定產(chǎn)品的纖維增強物件中的纖維取向分布。本發(fā)明的優(yōu)點是可以以提高的精度確定產(chǎn)品的纖維增強物件中的纖維取向分布。本發(fā)明的另一個優(yōu)點是可以更快地確定產(chǎn)品的纖維增強物件中的纖維取向分布。
雖然為了說明的目的具體描述了本發(fā)明，應(yīng)該理解這些細節(jié)只是為了說明的目的，本領(lǐng)域技術(shù)人員可以在其中做各種變型而不脫離本發(fā)明的范圍。
權(quán)利要求中使用的術(shù)語“包含”并不排除其它單元或步驟。單稱并不排除復(fù)指。
參考文獻
S.G.Advani(Ed.)Flow and Rheology in Polymer CompositesManufacturing，Elsevier，阿姆斯特丹(1994)
G.B.JefferyThe motion of ellipsoidal particles immersed in a viscousfluid，Proc.R.Soc.A 102，161-179(1922)
M.Junk和R.IllnerA new derivation of Jeffery′s equation，J.Math.Fluid Mech.8，1-34(2006)
F.Folgar和C.L.Tucker IIIOrientation behaviour of fibers inconcentrated suspensions，J.Reinf.Plast.Compos.3，98-119(1984)
C.L.Tucker和S.G.AdvaniProcessing of short-fiber systems，ch.6，pp.147in S.G.Advani(Ed.)Flow and Rheology in Polymer CompositesManufacturing(見以上參考文獻[1])
S.G.Advani和C.L.Tucker IIIThe use of tensors to describe andpredict fiber orientation in short fiber composites，J.Rheol.，751-784(1987)
S.G.Advani和C.L.Tucker IIIClosure approximations forthree-dimensional structure tensors，J.Rheol.，367-386(1990)
J.S.Cintra和C.L.Tucker IIIOrthotropic closure approximationsfor flow-induced fiber orientation，J.Rheol.39，1095-1122(1995)
E.Hairer、C.Lubich和G.WannerGeometric numerical integration，Springer，柏林(2002)
J.Linn、J.Steinbach、A.ReinhardtCalculation of the 3D fiberorientation in the simulation of the injection molding process forshort-fiber reinforced thermoplasts，ECMI 2000 Conference，Palermo(2000)
J.LinnExploring the phase space of the Folgar-Tucker equation，SIAM-EMS Conference，柏林(2001)
J.LinnOn the frame-invariant description of the phase space of theFolgar-Tucker equation，p.327-332in A.Buikis，R.
A.D.Fitt(Eds.)Progress in Industrial Mathematics at ECMI 2002，Springer(2004)
S.HessFokker-Planck equation approach to flow alignment in liquidcrystals，Z.Naturforsch.31A，1034ff.(1976)
M.DoiMolecular dynamics and rheological properties ofconcentrated solutions of rodlike polymers in isotropic and liquidcristalline phases，J.Polym.Sci.，Polym.Phys.Ed.19，229-243(1981)
M.Grosso、P.L.Maffetone、F.DupretA closure approximation fornematic liquid crystals based on the canonical distribution subspacetheory，Rheol.Acta39，301-310(2000)
M.
Simple models for complex nonequilibrium fluids，Phys.Rep.390，453-551(2004)
J.LinnThe Folgar-Tucker Model as a Differential Algebraic Systemfor Fiber Orientation Calculation，pp.215-224 in Y.Wang，K.HutterTrends in Applications of Mathematics to Mechanics，proceedings ofthe STAMM 2004 conference in Seeheim(德國)，Shaker(2005)
U.StrautinsInvestigation of fiber orientation dynamics within theFolgar-Tucker model with hybrid closure，碩士論文，數(shù)學(xué)系，凱澤斯勞滕大學(xué)(2004)
J.LinnDreidimensionale Vorausberechnung der anisotropenMaterialeigenschaften in thermoplastischen Spritzgusserzeugnissen(AnIso-3D)，Projektbericht für die MAGMA GmbH，Teil IIa(2002)
J.LinnDreidimensionale Vorausberechnung der anisotropenMaterialeigenschaften in thermoplastischen Spritzgusserzeugnisse(AnIso-3D)，Projektbericht für die MAGMA GmbH，Teil IIb(2002)
B.E.Ver Weyst、C.L.Tucker、P.H.Foss、J.F.O′GaraFiberorientation in 3D injection moulded featuresprediction andexperiment，Internat.Polymer Processing 14，409-420(1999)；
B.E.VerWeystNumerical predictions of flow-induced fiberorientation in three-dimensional geometries，博士論文，伊利諾伊大學(xué)厄本那-香檳分校(1998)
G.I.MarchukSplitting and Alternating Direction Methods，pp.197-462in P.G.Ciaret&J.L Lions(Eds.)Handbook of NumericalAnalysis，卷I，北荷蘭，Elsevier(1990)
K.W.MortonNumerical solution of convection-diffusion problems，Chapman&Hall，倫敦(1996)
R.J.LeVeque，Numerical Methods for Conservation Laws，
(1992)
G.Strang″On the construction and comparison of differenceschemes″，SIAM Journ.Num.Anal.5，506-517(1968)
M.G.Crandall和A.Majda″The method of fractional steps forconservation laws″，Math.Comp.34，285-314(1980)
H.V.Kojouharov、B.M.Chen″Nonstandard methods for theconvective transport equation with nonlinear reactions″，Numer.Meth.Partial Diff.Eq.14，467-485(1998)；″Nonstandard methods for theconvective-dispersive transport equation with nonlinear reactions″inR.E.Mickens(ed.)Applications of non-standard finite differenceschemes，minisymposium on nonstandard finite difference schemestheory and applications，SIAM年會，亞特蘭大GA，USA 1999，由新加坡World Scientific(2000)發(fā)表
H.Wang、X.Shi和R.E.Ewing″An ELLEM scheme formultidimensional advection-reaction equations and its optimal-ordererror estimate″，SIAM J.Numer.Anal.38，1846-1885(2001)
P.J.van der Houwen″Note on the time integration of 3Dadvection-reaction equations″，J.Comput.Appl.Math.116，275-278(2000)
W.Hunsdorfer，J.G.Venver″A note on splitting errors foradvection-reaction equations″，Appl.Numer. Math.18，191-199(1995)
S.V.PatankarNumerical heat transfer andfluidflow，HemispherePubl.Corp.(1980)
C.A.J.FletcherComputational techniques for Fluid Dynamics，卷IFundamental and General Techniques(第二版)，Springer(1991)
L.F.Shampine″Conservation laws and the numerical solution ofODEs″，Comp.Math.Applic.12B，1287-1296(1986)
L.F.Shampine″Linear conservation laws for ODEs″，Comp.Math.Applic.35，45-53(1998)
L.F.Shampine″Conservation laws and the numerical solution ofODEs，part II″，Comp.Math.Applic.38，61-72(1999)
J.Linn″Entwicklung eines Software-Moduls zur Berechnung derFaserorientierung in der Spritzgieβsimulation mit SIGMASOFT″，technical Report for the MAGMA GmbH(2001)
W.H.Press、S.A.Teukolsky、W.T.Vetterling、B.P.FlanneryNunlerical Recipes in Fortran 77The Art of Scientific Computing(第二版)，劍橋大學(xué)出版社(1992)
J.Stoer、R.BulirschIntroduction to Numerical Analysis(第三版)，Springer(2002)
D.H.Chung、T.H.Kwon″An invariant-based optimal fittingclosure approximation for the numerical prediction of flow-inducedfibre orientation″，J.Rheol.46，169-194(2002)
F.Dupret，V.Verleye，B.Languilier″Numerical prediction ofmoulding of short-fibre composite parts″，Proc.1st ESAFORM Conf.，291-294(1998)
G.H.Golub，H.A.van der Vorst″Eigenvalue Computation in the20th Century″，J.Comp.Appl.Math.123，35-65(2000)
權(quán)利要求
1.一種用于計算非球形粒子在宏觀水平上的取向統(tǒng)計的計算機實現(xiàn)的方法，該方法通過使用用于模擬注塑成型過程的模擬模型來進行所述計算，在所述注塑成型過程中，形成模擬區(qū)域的至少部分幾何形狀的模具型腔被懸浮液填充或部分地填充，所述懸浮液由包含大量非球形粒子的溶劑形成，其中，提供所述模擬區(qū)域的幾何形狀的數(shù)字表示或計算機模型，并且其中通過對至少部分所述模擬區(qū)域進行細分或者離散化而形成具有多個計算單元的網(wǎng)格，所述方法包括步驟
a)指定邊界條件；
b)設(shè)置初始條件；
c)求解針對所述模擬區(qū)域的至少部分單元的質(zhì)量、動量和能量的平衡方程，以獲得宏觀水平上的流體流、熱流和質(zhì)量傳遞；和
d)至少部分地基于所求解的平衡方程的結(jié)果求解非球形粒子取向動力學(xué)方程，以由此確定作為空間和時間的函數(shù)的宏觀水平上非球形粒子取向的改變。
2.根據(jù)權(quán)利要求1所述的方法，其中對于步驟d)，僅針對包含懸浮液的計算單元求解所述粒子取向方程。
3.根據(jù)權(quán)利要求1或2所述的方法，其中步驟c)還包括(cc)至少部分基于所求解的平衡方程的結(jié)果，確定所述流體或懸浮液的經(jīng)更新的自由表面或流前沿，其中所述自由表面將所述型腔的被所述懸浮液填充的單元與空單元分開。
4.根據(jù)權(quán)利要求3所述的方法，其中步驟(cc)還包括根據(jù)所述經(jīng)更新的流前沿更新所述邊界條件。
5.根據(jù)權(quán)利要求3或4所述的方法，還包括步驟
e)通過確定所述模具型腔是否被所述懸浮液的模擬注入填滿來確定所模擬的注塑成型過程是否結(jié)束；和
f)重復(fù)步驟c)、cc)和d)，直至所模擬的注塑成型過程結(jié)束為止。
6.一種用于在過程的模擬中描述非球形粒子的統(tǒng)計分布取向的方法，在所述過程中，用由包含大量非球形粒子的溶劑形成的懸浮液填充模具型腔，所述方法包括步驟
提供定義所述型腔的幾何形狀的三維計算機模型；
指定邊界條件；
基于所述模型對解域進行離散化以形成具有多個單元的網(wǎng)格；
針對所述解域的至少部分求解能量和流方程；
計算作為時間的函數(shù)的各個單元中的流條件和溫度條件；
計算非球形粒子取向的改變；和
描述作為時間的函數(shù)的所述各個單元中所述非球形粒子的取向的統(tǒng)計分布。
7.根據(jù)權(quán)利要求1至6中任意權(quán)利要求所述的方法，所述方法使用將所述流體流方程從所述粒子取向方程部分解耦的原理，由此獨立于所述粒子取向方程的計算進行所述流體流方程的計算，并且其中根據(jù)所求解的平衡方程的所述結(jié)果獲得局部流速和所述局部流速的空間梯度，所述局部流速和所述局部流速的梯度用作所述粒子取向方程的輸入。
8.根據(jù)權(quán)利要求7所述的方法，其中通過模擬模型對所述懸浮液的流進行建模，所述模擬模型通過在所有計算單元中、在離散時間瞬時tn處求解針對質(zhì)量、動量和能量的離散化輸運方程來產(chǎn)生所述懸浮液的狀態(tài)變量，包括流速U(r，t)及其空間梯度，所述計算單元位于空間中的離散點rm周圍，所述離散點包含在所述解域的被填充部分中。
9.根據(jù)權(quán)利要求1至8中任意權(quán)利要求所述的方法，其中所述模具型腔是注塑成型機的部分，并且所述成型機具有用于將所述懸浮液送入到所述模具型腔中的澆注系統(tǒng)和澆口，并且其中所述澆注系統(tǒng)和所述澆口的幾何形狀也形成所述模擬區(qū)域的部分。
10.根據(jù)權(quán)利要求1至9中任意權(quán)利要求所述的方法，其中所述非球形粒子是纖維和/或小片。
11.根據(jù)權(quán)利要求1至10中任意權(quán)利要求所述的方法，其中當(dāng)所述溶劑的粘度足夠大時，忽略慣性對所述非球形粒子的運動的影響。
12.根據(jù)權(quán)利要求11所述的方法，其中假設(shè)所述非球形粒子足夠小，從而能夠在粒子附近假設(shè)恒定的速度梯度。
13.根據(jù)權(quán)利要求1至12中任意權(quán)利要求所述的方法，其中利用分布函數(shù)在宏觀水平上統(tǒng)計地描述所述懸浮液中的粒子或纖維取向(FO)。
14.根據(jù)權(quán)利要求13所述的方法，其中由遞增階的粒子或纖維取向(FO)張量近似所述分布函數(shù)，并且其中在所述分布函數(shù)的近似中使用二階和四階粒子或纖維取向張量。
15.根據(jù)權(quán)利要求14所述的方法，其中在所述分布函數(shù)的近似中僅僅使用二階和四階粒子或纖維取向張量。
16.根據(jù)權(quán)利要求14或15所述的方法，其中近似解用于四階取向張量，通過使用閉合近似根據(jù)二階取向張量計算所述近似四階取向張量。
17.根據(jù)權(quán)利要求16所述的方法，其中所述閉合近似是混合閉合近似。
18.根據(jù)權(quán)利要求17所述的方法，其中所述混合閉合近似已經(jīng)被穩(wěn)定化以實現(xiàn)數(shù)值穩(wěn)定性。
19.根據(jù)權(quán)利要求18所述的方法，其中所述混合閉合近似的穩(wěn)定化包括使用標(biāo)量取向因子fSC，將所述標(biāo)量取向因子限制到
的數(shù)值范圍。
20.根據(jù)權(quán)利要求19所述的方法，其中通過在所計算出的fSC的值低于0時將所述值設(shè)置為0、并且在所計算出的fSC的值高于1時將所述值設(shè)置為1來將所述標(biāo)量取向因子的值限制到數(shù)值范圍
。
21.根據(jù)權(quán)利要求1至20中任意權(quán)利要求所述的方法，其中由其作為空間坐標(biāo)和時間坐標(biāo)的函數(shù)的取向向量給出的粒子或纖維的瞬時取向基于用于單個纖維的確定性取向動力學(xué)的Jeffery模型，還包括用于隨機地對纖維-纖維相互作用進行建模的隨機項。
22.根據(jù)權(quán)利要求21所述的方法，其中，使用包括用于對纖維-纖維相互作用進行建模的旋轉(zhuǎn)擴散系數(shù)的所述粒子或纖維分布函數(shù)的Fokker-Planck方程來對所述粒子或纖維的瞬時取向動力學(xué)進行建模。
23.根據(jù)權(quán)利要求22所述的方法，其中遵從Folgar和Tucker的方法，并且將所述旋轉(zhuǎn)擴散系數(shù)設(shè)置為與所述局部速度梯度的有效剪切速率成比例。
24.根據(jù)權(quán)利要求22或23以及權(quán)利要求14至21中任意權(quán)利要求所述的方法，其中，結(jié)合所述Fokker-Planck方程使用所述分布函數(shù)的纖維取向(FO)張量近似，以產(chǎn)生名為Folgar-Tucker方程(FTE)的、針對遞增階的FO張量的相互耦合方程的族，并且其中所要求解的粒子取向方程至少部分基于所述Folgar-Tucker方程或經(jīng)修改的Folgar-Tucker方程。
25.根據(jù)權(quán)利要求1至24中任意權(quán)利要求所述的方法，其中所要求解的粒子取向方程至少部分基于所謂的Folgar-Tucker方程(FTE)或經(jīng)修改的Folgar-Tucker方程。
26.根據(jù)權(quán)利要求1至25中任意權(quán)利要求所述的方法，其中利用所述Folgar-Tucker方程計算粒子取向和輸運的改變。
27.根據(jù)權(quán)利要求24、25或26所述的方法，其中將閉合近似應(yīng)用于所述FTE的二階FO張量方程以計算所述四階FO張量的近似解。
28.根據(jù)權(quán)利要求24至27中任意權(quán)利要求所述的方法，其中利用二階矩描述所述懸浮液中的非球形粒子的取向的統(tǒng)計分布。
29.根據(jù)權(quán)利要求28所述的方法，其中通過所述二階矩描述所述非球形粒子的取向的統(tǒng)計分布，所述二階矩是取向分布密度函數(shù)ψ(p)的二階取向張量。
30.根據(jù)權(quán)利要求27至29中任意權(quán)利要求所述的方法，其中將穩(wěn)定化混合閉合近似應(yīng)用于所述二階FO張量方程以計算所述4階FO張量的近似解，或者其中使用穩(wěn)定化混合閉合近似將所述四階取向張量作為所述二階張量的函數(shù)進行計算。
31.根據(jù)權(quán)利要求30所述的方法，其中所述混合閉合近似的穩(wěn)定化包括使用標(biāo)量取向因子，將所述標(biāo)量取向因子限制到
的數(shù)值范圍。
32.根據(jù)權(quán)利要求30或31所述的方法，其中通過經(jīng)由添加懲罰項修改所述Folgar-Tucker方程來動態(tài)地穩(wěn)定所述二階張量的跡，即所述二階張量的對角元的和
33.根據(jù)權(quán)利要求32所述的方法，其中特別地選擇所述懲罰項的函數(shù)形式，使得將所述二階取向張量的跡Tr[aij(2)]近似地保持在其要求值1。
34.根據(jù)權(quán)利要求33所述的方法，其中通過選擇所述懲罰項來將所述取向張量近似地保持在其要求值1，以使得具有單位跡的所有對稱3×3矩陣的集合變?yōu)樗鼋?jīng)修改的Folgar-Tucker方程的穩(wěn)定積分流形。
35.根據(jù)權(quán)利要求30至34中任意權(quán)利要求所述的方法，其中通過將所述標(biāo)量取向因子的值歸到其理論上所容許的0和1之間的范圍來穩(wěn)定所述混合閉合近似，在利用所述混合閉合近似的所述四階FO張量的計算中使用所述標(biāo)量取向因子。
36.根據(jù)權(quán)利要求31或35所述的方法，其中通過在所計算出的fSC的值低于0時將所述值設(shè)置為0、并且在所計算出的fSC的值高于1時將所述值設(shè)置為1來將所述標(biāo)量取向因子fSC的值限制到所述數(shù)值范圍
。
37.根據(jù)權(quán)利要求24至36中任意權(quán)利要求所述的方法，其中用于所述纖維取向的計算的時間積分方法局部地適應(yīng)于局部剪切速率水平。
38.根據(jù)權(quán)利要求37所述的方法，其中所述適應(yīng)包括
(a)時間步長隨剪切速率的增大而減??；和/或
(b)所述時間積分方法的階隨剪切速率的增大而提高；和/或
(c)隨著剪切速率的增大，將所述時間步局部地分成若干個更小的時間步。
39.根據(jù)權(quán)利要求38所述的方法，其中以這樣的方式使用所述適應(yīng)(a)、(b)和(c)中的至少一個由流計算給出的全局時間步上的誤差對于所有被填充的單元近似相同，其不依賴于所述剪切速率的局部變化，即使在所述局部變化非常大的情況下也是如此。
40.根據(jù)權(quán)利要求24至39中任意權(quán)利要求所述的方法，其中利用不變量來控制如下限制由所述纖維取向矩陣的第二不變量和第三不變量的數(shù)值將所述Folgar-Tucker方程的數(shù)值解限制到其允許的域或相空間。
41.根據(jù)權(quán)利要求1至40中任意權(quán)利要求所述的方法，其中所述懸浮液的壓縮性效應(yīng)包括在所述模型中，即使在所述壓縮性效應(yīng)非常小的情況下也是如此。
42.根據(jù)權(quán)利要求24至41中任意權(quán)利要求所述的方法，其中將所述Folgar-Tucker方程分成純粹取向輸運分量和旋轉(zhuǎn)取向動力學(xué)分量所述Folgar-Tucker方程是形式的反應(yīng)對流類型的雙曲型PDE的非線性耦合系統(tǒng)。
43.根據(jù)權(quán)利要求42所述的方法，其中通過以交替的方式使用根據(jù)算子分裂過程的兩個方程來執(zhí)行所述數(shù)值積分。
44.根據(jù)權(quán)利要求43所述的方法，其中通過對稱算子分裂或通過簡單算子分裂來進行所述算子分裂。
45.根據(jù)權(quán)利要求1至44中任意權(quán)利要求所述的方法，其中，不管所述非球形粒子在真實的模具填充過程中將具有的分布，在所述模擬開始時所述非球形粒子的初始取向是隨機化分布。
46.根據(jù)權(quán)利要求1至45中任意權(quán)利要求所述的方法，其中在所述纖維取向的初始化時忽略在所述自由表面處的粒子的噴泉流效應(yīng)。
47.根據(jù)權(quán)利要求1至46中任意權(quán)利要求所述的方法，其中在所述填充模擬的每個時間步時根據(jù)來自相鄰單元的質(zhì)量輸運、通過加權(quán)平均來初始化所述自由表面處新填充的單元中的纖維的纖維取向分布。
48.根據(jù)權(quán)利要求24至48中任意權(quán)利要求所述的方法，其中使用簡化符號(CN)以向量形式重新用公式表示所述Folgar-Tucker方程(16)。
49.根據(jù)權(quán)利要求24至48中任意權(quán)利要求所述的方法，其中利用所述Folgar-Tucker方程(16)的右側(cè)的代數(shù)結(jié)構(gòu)以很少的代數(shù)運算實現(xiàn)穩(wěn)定化混合閉合。
50.根據(jù)權(quán)利要求49所述的方法，其中通過識別所述Folgar-Tucker方程(16)的右側(cè)中的共同子表達式來使確定所述纖維取向所需要的代數(shù)運算的數(shù)量最小化。
51.根據(jù)權(quán)利要求49或50所述的方法，其中通過將所述混合閉合的定義公式代入到右側(cè)函數(shù)中來重新用公式表示所述Folgar-Tucker方程(16)，用于對所得到的項進行解析計算和隨后的重組，以實現(xiàn)對于高效估算最優(yōu)的項結(jié)構(gòu)。
52.根據(jù)權(quán)利要求48至51中任意權(quán)利要求所述的方法，其中在所述重新用公式表示的Folgar-Tucker方程(16)中包括控制項用于動態(tài)跡穩(wěn)定化。
53.根據(jù)權(quán)利要求48至52中任意權(quán)利要求所述的方法，其中使用CN以最小數(shù)量的運算來估算所述重新用公式表示的右側(cè)函數(shù)用于最優(yōu)效率。
54.根據(jù)權(quán)利要求48至53中任意權(quán)利要求所述的方法，還包括檢查所述FTE的數(shù)值解是否在所述第二不變量和第三不變量的控制的可容許的域內(nèi)。
55.根據(jù)權(quán)利要求48至54中任意權(quán)利要求所述的方法，還包括使用“跡尺度改變”以解釋通過DTS而僅具有近似單位跡的FO矩陣。
56.根據(jù)權(quán)利要求48至55中任意權(quán)利要求所述的方法，包括相空間投影的算法描述。
57.根據(jù)權(quán)利要求56所述的方法，還包括相空間投影和其它所提出的用于所述FTE的數(shù)值積分的技術(shù)的結(jié)合以將所述數(shù)值解限制到所述可容許的域。
58.根據(jù)權(quán)利要求27至57中任意權(quán)利要求所述的方法，其中將作為所述Folgar-Tucker方程(FTE)或者經(jīng)修改的Folgar-Tucker方程的解的解軌跡或FO張量的集合限制到粒子或者纖維取向張量或矩陣的相空間集合MFT，其中矩陣的所述相空間MFT集合是半正定且具有單位跡的實對稱3×3矩陣的向量空間的有界凸子集，由此可以將所述FTE看作是差分代數(shù)系統(tǒng)(DAS)。
59.根據(jù)權(quán)利要求58所述的方法，其中所述FTE的求解包括監(jiān)測步驟以確定FO張量或者矩陣的所提出的解是否是所允許的解，其中所述張量或者矩陣為半正定并且具有單位跡，并且如果所提出的解不是所允許的解則執(zhí)行相空間投影步驟，其中將所提出的解投影到所述相空間MFT上。
60.根據(jù)權(quán)利要求59所述的方法，其中通過將所提出的FO矩陣映射到所述相空間MFT中的最近的矩陣上來執(zhí)行所提出的FO解矩陣的投影。
61.根據(jù)權(quán)利要求58或59所述的方法，其中所述監(jiān)測步驟是不變量監(jiān)測步驟，其包括確定所述FO張量或者矩陣的第二不變量和第三不變量的數(shù)值。
62.根據(jù)權(quán)利要求60或61所述的方法，其中將所提出的FO矩陣映射到所述相空間MFT中的最近的矩陣上包括使用所提出的FO矩陣的特征值的跡尺度改變。
63.根據(jù)權(quán)利要求7以及權(quán)利要求24至62中任意權(quán)利要求所述的方法，其中求解所述FTE包括使用對稱算子分裂，其中所述對稱算子分裂為對稱包括方法，其包括擴展步驟、預(yù)旋轉(zhuǎn)步驟、對流步驟和后旋轉(zhuǎn)步驟。
64.根據(jù)權(quán)利要求24至63中任意權(quán)利要求所述的方法，其中求解所述FTE的方程包括使用所述FO張量或矩陣的所有六個元素。
65.根據(jù)權(quán)利要求24至64中任意權(quán)利要求所述的方法，其中在時間間隔內(nèi)求解所述FTE的方程包括使用經(jīng)縮放的流速梯度和經(jīng)縮放的時間坐標(biāo)。
66.根據(jù)權(quán)利要求34至65中任意權(quán)利要求所述的方法，其中求解所述FTE的方程包括使用結(jié)合了通過控制項進行的動態(tài)跡穩(wěn)定化的流受控時間積分。
67.根據(jù)權(quán)利要求1至66中任意權(quán)利要求所述的方法，其中所述模擬的結(jié)果用于改進要通過被模擬的所述注塑成型過程制造的物件的特性。
68.根據(jù)權(quán)利要求67所述的方法，其中所述模擬的結(jié)果用于改進所要制造的物件的強度、剛性或穩(wěn)定性特性。
69.根據(jù)權(quán)利要求67所述的方法，其中所述模擬的結(jié)果用于提高形狀穩(wěn)定性，以降低所要制造的物件的變形量。
70.根據(jù)權(quán)利要求1至69中任意權(quán)利要求所述的方法，還包括通過所述粒子取向的統(tǒng)計分布來表征流計算的宏觀水平上的局部纖維取向狀態(tài)。
71.一種計算機可讀介質(zhì)上的計算機產(chǎn)品，其包括用于執(zhí)行根據(jù)權(quán)利要求1至70中任意權(quán)利要求所述的方法的軟件算法。
72.一種用于最佳的鑄造產(chǎn)品的系統(tǒng)，其包括產(chǎn)品鑄造機械和與其連接的計算機實現(xiàn)的優(yōu)化系統(tǒng)，所述優(yōu)化系統(tǒng)包括根據(jù)權(quán)利要求71所述的計算機產(chǎn)品。
全文摘要
一種用于在過程的模擬中描述非球形粒子的統(tǒng)計取向分布的方法和裝置，在該過程中用包含大量非球形粒子的懸浮液填充模具型腔?？梢詫⑺龇椒ê脱b置應(yīng)用于制造纖維增強模制聚合物部件的注塑成型過程的分析或制造纖維增強金屬產(chǎn)品的金屬鑄造過程的分析。這些分析的結(jié)果可以用于確定部件的張力和翹曲的方面，并且用于優(yōu)化在制造過程中使用的過程條件。
文檔編號G06F17/50GK101754845SQ200880016110
公開日2010年6月23日 申請日期2008年7月1日 優(yōu)先權(quán)日2007年7月2日
發(fā)明者約阿希姆·林, 馬蒂亞斯·莫格 申請人:馬格馬鑄造工藝有限公司