歐拉法的案例
流固耦合（FESIM有限元分析）
流固耦合包含一般耦合、任意偶合，且采用拉格朗日法與歐拉法分別描述固體與流體的運動。拉格朗日的元素節點依附在材料上，節點隨著材料質點作運動，故各物理量也作用在節點上隨材料流動而變化。相反，除任意耦合外，歐拉元素網格與節點不隨時間而變，其物理量雖也作用在歐拉元素節點上，但對于通過歐拉元素面的各時間的質量、動量與能量的進與出，加之模擬，即模擬元素面的材料流，而不模擬各材料質點的時間歷程。因為對一般固體材料，要模擬各材料質點的時間歷程，因此大多用拉格朗日法。而對于流體不需要模擬材料質點的時間歷程，故采用歐拉法，歐拉法需用三維的計算域、三維的體元素與通用材料。此外，歐拉法容許一個元素內含有兩種以上的材料，這就是模擬計算材料流的擴散與混合行為。
拉格朗日法與歐拉法是對運動現象的兩類不同的數學描述，可說是分別對材料質點流與空間流之描述。拉格朗日法與歐拉法之元素網格可在同一計算模型內，但拉格朗日法的元素與歐拉法的元素分別擁有節點，只采用介面(interface)，稱為耦合面，才能將兩者連結在一起；否則，縱使兩者在空間內相互重疊，也彼此不相干，即忽視對方之存在。
1 拉格朗日法
對固定的坐標系而言，拉格朗日元素的節點可相對地運動。因節點系附在材料上，故材料連續體之節點系一起隨著材料質點流而運動。各拉格朗日元素的質量是不變量(invariant)，但其元素體積可隨時間而改變。此外，速度、壓力強度或質量密度等物理量系作用在拉格朗日元素的節點上，因此，各物理量系隨著材料流(material flow)而改變。因對固體材料之行為， 較須追蹤各材料質點之時間歷程，故適宜采用拉格朗日法。拉格朗日法也適宜用以分析材料破壞(failure)或應變硬化(strain hardening)問題。
展開 耦合歐拉-拉格朗日（CEL）法攪拌摩擦焊接模擬
20170929205500.gif
采用耦合歐拉-拉格朗日法對攪拌摩擦焊接攪拌頭下扎過程進行Abaqus數值模擬。
CFD理論|基本方程（1）
（1）拉格朗日法
該方法著眼于流體質點，把流體質點的物理量表示為拉格朗日坐標與時間的函數。拉格朗日法跟蹤的是流體質點，因此其坐標（a,b,c）不隨時間(T)的變化。若以 表示流體質點的某一物理量，其拉格朗日法的數學表示方式為：
(a,b,c)更像是對流體質量的標號，如果t時刻的質點的位置以r表示，則：
表示拉格朗日坐標為(a,b,c)的流體質點在t時刻處于r，即空間點(x,y,z)的位置。
（2）歐拉法
歐拉法著眼于空間點，也叫空間描述法，即流體的物理量隨空間點及時間變化。把流體物理量表示為歐拉坐標及時間的函數。當某時刻物理量在空間分布確定，我們就是物理量在此空間形成了一個場，也就是說歐拉法實際上描述了物理量的場。歐拉的數學表達式為：
(x,y,z)就是空間的坐標。例如流體速度可以表示為：
它表示空間點(x,y,z)上t時刻的流體速度。
（3）隨體導數
拉格朗日法描述中，物理量隨體導數就是質量(a,b,c)的物理量對時間的導數，比如說速度及加速度的表達：
在歐拉法描述中， 并不表示隨體導數，它只是表示物理量在空間點（x,y,z)隨時間的變化規律，因為隨體導數針對的對象是流體質點，而不是空間點，因此還需要考慮空間點(x,y,z)隨時間的變化。所以物理量隨體導數的歐拉表示方法為：
簡化后可以表示為：
這是流體力學十分重要的基本公式。方程右邊第一項為遷移加速度/局部導數，表示該空間點物理量隨時間的變化規律；第二項為當地加速度/位變導數項，它是物理量不定常性引起，由空間位置變化引起的物理量變化律。
雷諾輸運定理
隨體導數后，我們如何去描述流體性質的變化，這就需要用到雷諾輸運定理。在介紹這個定理前，我們首先要知道控制體與系統。
展開 ABAQUS UMAT for cohesive element shear behaviour ￥150
本文利用 cohesive element 模擬兩個粘接的物體受剪破壞，利用的本構模型是Coulomb模型，有限元實現算法是帶有返回功能的前歐拉法（Crisfield在1991出版的Non-linear finite element analysis of solid and structure中的第6章有提到），軟件是abaqus，子程序為UMAT。
這種有限元模型主要應用在在膠合破壞的預測，比如磚結構的水泥砂漿、纖維復合材料加固結構中復合材料的剝落、多層玻璃的脫膠以及夾層板的滑移等等。可以說，與膠有關的結構都可以試著利用cohesive單元模擬，所以開發cohesive單元的本構模型是很多領域的熱點，因為相比起固體單元，cohesive單元的發展相對較晚。
和tension cut-off那篇一樣，我用前歐拉法推導出的tangent stiffness和Lourenco（1997）用后歐拉法推導出的tangent stiffness是一樣的。
Coulomb的返回算法：
視頻結果詳見：https://zhuanlan.zhihu.com/p/113150354
模擬結果：
應力：
應變：
步驟1：part （2D 模型）
步驟2：material
步驟3：assembly
步驟4：step
步驟5：interaction
步驟6：load
步驟7：mesh
步驟8：結果
展開 
ABAQUS UMAT for cohesive element tension behaviour ￥150
本文利用 cohesive element 模擬兩個粘接的物體受拉破壞，利用的本構模型是tension-cut off模型，有限元實現算法是帶有返回功能的前歐拉法（Crisfield在1991出版的Non-linear finite element analysis of solid and structure中的第6章有提到），軟件是abaqus，子程序為UMAT。
這種有限元模型主要應用在在膠合破壞的預測，比如磚結構的水泥砂漿、纖維復合材料加固結構中復合材料的剝落、多層玻璃的脫膠以及夾層板的滑移等等。可以說，與膠有關的結構都可以試著利用cohesive單元模擬，所以開發cohesive單元的本構模型是很多領域的熱點，因為相比起固體單元，cohesive單元的發展相對較晚。
這里有意思的是，我用前歐拉法推導出的tangent stiffness和Lourenco（1997）用后歐拉法推導出的tangent stiffness是一樣的。
tensile cut-off的返回算法：
模擬的視頻結果詳見：
https://zhuanlan.zhihu.com/p/113143055
模擬結果：
mises應力圖
ABAQUS建模：這是一個二維的單元，兩個彈性材料物體被一層cohesive單元（厚度為10）連接著。
展開 十六、DPM模型-顆粒流動
1 概念介紹
首先我們介紹一下拉格朗日法和歐拉法，理解起來很簡單，拉格朗日法是以某一質點的運動作為研究對象，觀察這一質點在流場中由一點移動到另一點時，其運動參數的變化規律；歐拉法以某一流場區域作為研究對象，研究各時刻質點在流場中的變化規律。
顯然，拉格朗日法更適用于描述顆粒運動，而歐拉法更適用于描述流體運動。DPM模型就是基于這兩種方法進行流體相和顆粒相的模擬，它使用歐拉法描述流體運動，使用拉格朗日法描述顆粒運動。
DPM適用條件：DPM模型只適用于顆粒相體積分數小于10%，同時不考慮顆粒體積。不考慮顆粒和顆粒之間的相互作用力，但可以考慮顆粒和流體之間的相互作用。
2 模型描述
本例的模型采用三通管模型，如圖所示。模型有兩個入口和一個出口，分別為INLET_Y、INLET_Z和OUTLET，含顆粒物的空氣從INLET_Z進口流入計算域內，最后經OUTLET流出。
3 導入網格
使用Workbench打開工程文件，文件在本文末尾鏈接資源內。
4 Scale網格尺寸
Scale修改網格尺寸。如圖所示。
確保計算域尺寸是我們所需要的。本例中x方向尺寸-0.038~0.038m，y方
向-0.2~0.2m，z方向-0.038~0.2m
5 設置求解器
選擇壓力基(pressure-based)求解器，同時選擇穩態模擬，不考慮考慮重力。
6 設置計算模型
多相流模型不必打開，不考慮能量方程。湍流模型為標準的k-e模型，增強的壁面函數，打開Discrete Phase模型。
不勾選Interaction with Continuous Phase，不考慮顆粒相和流體相之間的相互作用。其余選擇均保持默認。
展開 你是學流體力學的？去，給我倒杯水
歐拉法 VS 拉格朗日法
說完了實質，再看看形式——流體力學的兩種主要描述方式：拉格朗日法和歐拉法。
在經典場論中，拉格朗日法（又稱體系法）是研究流場內個別流體質點在不同時刻的位置、流速、壓力等參數的變化。也就是用隨時間的變化來描述流體質點的運動參數，各質點運動狀況的總和就構成了整個流體的狀態。
歐拉法（又稱控制體法）是研究整個流場內不同時刻，不同位置上的流體質點的參數。它把注意力集中在選定的空間點上而不是選定的流體質點上。因此，歐拉法確定的是占據流場中的全部流體質點的瞬時流動參數。
為了更形象的對比歐拉法和拉格朗日法，小編整理了上面這幅圖：拉格朗日法就好像是跟在鴨子后面劃船，而歐拉法則更像是在站在橋上數多少鴨子游過去。需要特別說明的是，最早提出這兩類研究流體方法的人都是歐拉。
N-S方程 VS 玻爾茲曼方程
前面鋪墊了這么多，目的還是把玻爾茲曼這尊大神娓娓的請出來。當然在此之前，還有必要先來捋一捋經典力學描述流體的方式。
展開 在流體力學里，有兩種描述流體運動的方法：歐拉（Euler)和拉格朗日（Lagrange)方法。歐拉法描述的是任何時刻流場中各種變量的分布，而拉格朗日法卻是去追蹤每個粒子從某一時刻起的運動軌跡。
在一個風和日麗的午后，YC坐在河岸邊看河水流，恩，她總是很閑。如果YC的位置不動，她在自己目光能及的河面上劃出一塊區域，數某一時刻經過的船只數，如果可能的話，再數數經過的魚兒數；當然，如果手頭有些儀器，她可以干干正事，比如測測水流的速度、水的壓力、水的溫度等，由此得到每一時刻這一河流區域水流各物理量的分布。那么YC是在用歐拉方法研究流體。
這時，YC忽然看到一條船上坐著她的初戀情人，雖然根據陳安對初戀情人的定義，YC根本沒有初戀情人。現在假設她有，天哪，他們有20年沒見面了，他還欠她20元呢，不能放了他。于是YC記下第一眼看到初戀情人的時間，并迅速測出此時船的位置和速度，然后撒腿追去。假設這條船是順水而下，船的速度即是水流的速度。每隔一個時間點，她便測一下船的速度和位置。為了曾經的愛情，還有那不計利息的20元，她越過山崗，淌過小溪，直到那條船離開了她的視線。于是，她得到了這條船在河流中的運動軌跡。YC此時所用的研究方法就是拉格朗日法。
Understood?
而在一些復雜的兩相流動問題里，比如粒子在流場中運動的問題，我們關注的是粒子的運動軌跡，因此，我們可以用拉格朗日方法追蹤粒子在流場中的運動，同時，用歐拉方法來計算流場的各物理量。
在許多工程領域，都有纖維在流場中運動的問題。如果將纖維在流場中的運動視為兩相流動，必須為纖維作一些改變，因為它不同于一般的剛性粒子。它細長，細長到你無法用一個粒子來代表一根纖維；它柔，柔得自己的每一部分可以相對于其他部分發生變形。
展開 abaqus耦合的歐拉-拉格朗日單元技術建模分析法(含源文件) ￥5
1. 分別為Eulerian domain和Lagrangian domain建立兩個
建立Lagrangian domain的Part，類型設置為Discrete rigid，并設置Reference Point。建立Eulerian domain的Part，類型設置為Eulerian，要注意Eulerian domain和Lagrangian domain要保證有重疊的部分，這是一種弱耦合，數據在兩個區域間拋來拋去，所以網格要有重疊部分。
2. 定義水的材料屬性
選擇狀態方程模型EOS中Us-Up，設置聲速c0=1483m/s；密度為1000kg/m3；粘度為0.001kg/ms。并把截面屬性賦給Eulerian domain。
展開 ABAQUS UMAT后歐拉法實現cohesive單元的拉-剪混合本構 ￥1500
本文介紹了如何用后歐拉算法模擬2D cohesive單元的復雜破壞，采用的本構模型是hyperbolic曲線。這種曲線可以將受拉破壞和受剪破壞耦合在一起，很好地反應了膠合接觸面地特性。
主要參考文獻為：Caballero, A., Willam, K.J. and Carol, I., 2008. Consistent tangent formulation for 3D interface modeling of cracking/fracture in quasi-brittle materials. Computer Methods in Applied Mechanics and Engineering, 197(33-40), pp.2804-2822.
Willam, K.J 和 Carol, I. 長期致力于模擬巖石、混凝土和磚結構地開裂破壞。
本文的很多數學推導是通過軟件Wolfram Mathematica獲得的。
本文所用的算法是Caballero（2008）中的算法1：monolithic iteration strategy without substepping。
具體建模過程和建模結果見知乎文章：https://zhuanlan.zhihu.com/p/113538156
umat子程序和代碼對應的
詳細解答見附件。
展開 Abaqus金屬切削仿真
從計算力學的角度，切削仿真大致可以分為拉格朗日法、歐拉法、無網格法等。拉格朗日法是處理連續介質力學的經典有限元方法，現代的拉格朗日法也通過結合ALE、XFEM、單元剛度折減等技術手段來解決大變形、材料失效問題；歐拉法是計算流體力學的常用方法，也可用于固體力學中的大變形問題；無網格法包括SPH、DEM等，SPH方法也常用于切削過程中的大變形仿真。
ALE
CEL
SPH
高速銑削仿真
1、材料參數定義
通用參數：
材料密度
材料機械性能參數：
彈塑性階段（沒有損傷）
初始損傷準則（損傷起始）
損傷演化準則（剛度折減-材料失效）
材料熱力學性能參數：
導熱系數
線膨脹系數
比熱容
非彈性變形能耗散比
材料使用的是合金結構鋼20NiCrMo5，塑性階段我使用了J-C模型。
2、網格與單元
網格劃分要足夠細才會有切屑，刀具切削區域局部加密防止接觸穿透，單元類型選擇熱-位移耦合單元，要定義單元刪除和狀態輸出。
3、仿真結果
應力
溫度
等效塑性應變與切削力曲線
下期文章提供切削inp文件下載
鋁合金切削（未考慮溫度）
合金鋼切削（考慮溫度）
展開 
簡述極端變形問題的數值模擬
</span></p><p><span style="color: rgb(0, 0, 0);">由于背景網格本身不攜帶任何物質信息，且物質點和背景網格之間也沒有相對運動，這就避免了歐拉法中的對流項的計算，而且容易跟蹤物質界面，同時由于物質點攜帶了完整的物質信息，每次更新完狀態后就可以拋棄變形后的背景網格，避免了網格畸變問題的產生，實際的計算效率更高。</span></p><p class="ql-align-justify"><span style="color: rgb(0, 0, 0);">物質點方法綜合了</span><strong style="color: rgb(0, 0, 0);">拉格朗日法和歐拉法的優點，能夠很好的處理極端變形問題，例如沖擊侵徹、爆炸損傷、斷裂破碎、流固耦合等涉及結構和材料極端變形模擬</strong><span style="color: rgb(0, 0, 0);">；而對于小變形假設的靜力、擬靜力或低速動力學問題，相同插值階次的物質點法在計算精度上低于傳統有限元方法。因此，針對這類問題，采用有限元方法更合適。</span></p><p class="ql-align-justify"><strong style="color: rgb(0, 0, 0);">四、案例介紹</strong></p><p class="ql-align-justify"><span style="color: rgb(0, 0, 0);">物質點方法到底好不好用，相信大家通過上述理論知識已有初步了解。而該方法在實際應用中能否有效解決問題，我們不妨從以下幾個案例中尋找答案。
展開 fluent在氣溶膠傳播中的應用
使用數值模擬時，氣溶膠的模擬可以使用歐拉法和拉格朗日法，歐拉法將顆粒看做連續介質，在有限體積法基礎上構建濃度守恒方程，對氣溶膠顆粒的處理與對流相的處理方法類似；拉格朗日法將顆粒看做離散相，跟蹤每一個例子的運動軌跡。對例子的運動軌跡進行統計處理后，拉格朗日法也可以進行顆粒濃度及其他參數的計算。因顆粒濃度低，顆粒相和流體相可以看做單向作用，即流體主相影響顆粒相，而顆粒相對主相的影響忽略不計。
Fluent中提供了VOF、歐拉型和混合物型歐拉-歐拉多相流模型以及DPM歐拉拉格朗日模型，因氣溶膠在室內體積分數小于10%，我們使用DPM模型。
我們首先建立房間模型，尺寸為0.4X0.5X0.3m。右側窗戶為0.3X0.18m，下邊緣距離地面0.09m，4、6、5、7、8、10的尺寸均為0.04X0.01m，其中4、6是送風口；5、7、8、10是回風口。送風口9的尺寸為0.05X0.08m。各風口與房間的相對位置如圖；在房間中部我們使用兩個簡化的人體模型，他們距離0.1m，人體1對人體2持續講話，人體2暴露在飛沫氣溶膠顆粒中。
使用如下網格參數劃分出進出口及人體進行局部加密的網格。
為各房間入口、出口及人的嘴巴創建如下的邊界條件。
因為溫度對氣溶膠的傳播有影響，所以我們需要打開能量方程。在紊流模型中選擇標準k-e模型，根據進出口的尺寸，計算水力直徑為0.16m，打開離散相模型，設置離散相注入類型為surface，顆粒直徑為5μm，設置離散相濃度0.085μg/s、密度為1000kg/m。在邊界條件中設置人體1噴出的顆粒速度為3.9m/s，溫度35°，進風口速度均為2m/s，溫度31°。出風口均為自由出口，溫度35°。
展開 有限元理論基礎及Abaqus內部實現方式研究系列13：顯式和隱式的區別 ￥1
1.1 前向歐拉（Forward Euler）和后向歐拉（Backward Euler）
前向歐拉和后向歐拉分別是顯式和隱式的一個典型計算方法，本文將引用這兩個方法來盡量直觀地展現顯式求解和隱式求解的區別。
前向歐拉：
fn+1 (x) = fn (x) + h*f'n (x)
后向歐拉：
fn+1 (x) = fn (x) + h*f'n+1 (x)
1.2 算例
現在以一個微分方程算例為例：
y'(x) = -y+x+1
且有y(0) = 1。
1、前向歐拉法
使用前向歐拉法，可得：
yn+1 = yn + h*(-yn +xn +1)
顯然，這個式子不需要做任何額外運算就能從yn出yn+1，因此每步迭代計算量小。但h在一個最大限制，具體見下面的推導過程：
yn+1 = (1-h)yn + h*(xn +1)
可得
yn+1 = (1-h)2 yn-1 + h*(xn +1) + h*(xn-1 +1)
則得
yn+1 = (1-h)n +O(s2 )
當h<2時，y收斂，但h>=2時，y將發散。
2、后向歐拉法
yn+1 = yn + h*(-yn+1 +xn+1 +1)
這個式子不能直接得出yn+1, 必須做進一步計算得到
yn+1 = (yn + h*(xn+1 +1)) / (1+h)
當h>0時，y無條件收斂。
3、結果比較
假設n=30
（1）當h = 1.9時，計算結果如下圖所示,顯式和隱式都和理論接近：
（2）當h = 2.1時，計算結果如下圖所示，顯式求解發散了：
==總結==
本文簡單介紹了顯式和隱式的區別，本質在于采用不同的物理學平衡方程，因此在不同的物理學問題也有不同的表現。
展開 SPICE?電路仿真原理
這就意味著盡管電容事實上已經穩定了，仿真器仍然需要小的時間步長來保持它的電壓穩定，在數學術語里就是指此方法沒有剛性穩定，這種方法稱為前向歐拉積分法。SPICE 并沒有采用此方法，這里簡單介紹一下它可以幫助理解為何 SPICE 采用其他更復雜的方法來處理，其中最簡單的就是后向歐拉積分法。
我們知道 i0 ≈ C（v1-v0）/Δt， i1 ≈ C（v1-v0）/Δti1 并不是已知的，不能直接得到 v1 的表達式，意味著不能將電容。當作一個固定的電壓源來處理，而用一個電流源并聯一個電阻來替代電容：Δt 時 :i1=v1/req+ieq
i1 ≈ v1xC/Δt-v0*C/Δt 當 1/req=C/Δt，ieq= ﹣ v0*C/Δt 時，兩個方程等效。
因此我們可以用阻值為 Δt/C 的電阻并聯電流為﹣ v0*C/Δt 的電流源來替代電容。
這就是后向歐拉法，再對比之前提到的當電容電壓為 4.999V，time step 為1ms 時的情況：req=Δt/C=1e-3/1e-9=1e6。ieq= ﹣ v0*C/Δt= ﹣ 4.999x1e-9/1e-3= ﹣ 4.999e-6
因此計算的 v1 值為 4.999999001，精度很高，是剛性穩定的。SPICE 一般是在第一個時間點和之后的每個斷點時使用此方法處理，發生在波形的轉折點處，如脈沖的起始處。默認設置下，SPICE 通常采用另一種方法：梯形法則，也是剛性穩定的，同時具有更高的精度。
展開 