Newton-Raphson方法的案例
數學模型與非線性的定義——《非線性計算與多物理場耦合》系列課程之一
分別圍繞下面三個問題展開：
1.實際物理問題與數學模型之間的關系，怎么去建立或定義一個有效的數學模型，其與有限元方法的關系是什么？
2.我們為什么需要考慮非線性，非線性的數學關系式是什么，在有限元算法中體現在什么地方？
3.怎么運用基礎的Newton-Raphson方法去求解非線性方程組？
在視頻的中間穿插講述了本系列課程的基本框架，也就是一步一步非線性研究的每一個遞進關系的知識點，帶大家一步一步掌握非線性計算的相關知識。
此課附件包含兩個基于Julia寫的兩個代碼（Julia的安裝與基本操作視頻看完主頁的julia課程），PPT和完整視頻（免費完整視頻在我主頁課程里面），免費分享給大家，希望有興趣，覺得此視頻還有點用的同學關注我，后續會有更加精彩的內容。
Share1.zip
第一課 .pdf
展開 SIGMA/W 專業的應力變形有限元分析軟件
對非線性材料的土體模型，SIGMA/W軟件用Newton-Raphson方法來解迭代方程，土體參數在每一迭代步上都要更新，直至得到收斂解。
SIGMA/W軟件可與SEEP/W軟件相結合分析土體的完全固結問題。由SEEP/W軟件計算載荷作用下的瞬時孔隙水壓力，而由SIGMA/W軟件計算孔隙水壓力產生的變形。
用于土體結構內部相互作用的梁結構和桿單元。
回填或開挖時的分段載荷。
3.與其他應用軟件的結合：
1）SIGMA/W軟件計算出的應力可用于SLOPE/W軟件或QUAKE/W軟件中：
在SLOPE/W軟件中用有限元方法計算出的應力值，與由變形分析中得到的應力值一樣，用這些應力值就可以對一些嚴格的穩定性問題進行分析了。此外，在QUAKE/W軟件的地震動力學分析中，用戶可以將SIGMA/W軟件計算出的應力作為初始應力分布值。
2）SIGMA/W軟件計算出的孔隙水壓力可用于SLOPE/W軟件或QUAKE/W軟件中：
在SIGMA/W軟件中，在如回填等穩定載荷作用下產生的超孔隙水壓力可以代入SEEP/W軟件中研究地基中的超孔隙水壓力的消散所需時間。用戶可以用SLOPE/W軟件來分析建造過程這些附加應力對穩定性的影響，以便用戶決定分步加載的必要性。
展開 案例49-鋼筋混凝土板的載荷極限分析
• 使用初始Newton-Raphson非線性解方法更好地捕捉不穩定點。
• 與載荷控制分析相比，通過位移控制分析可以更容易地跟蹤剛度損失后的結構行為；然而，如果不穩定區域值得關注，并且需要進行載荷控制分析，則考慮使用弧長法（ARCLEN）。
使用弧長法，在大約610 kN的載荷和5.6 mm的撓度下確定了不穩定區域。結果與圖49.4所示的分析結果一致，驗證了計算的載荷極限。
參考文獻
Eurocode 2: Design of Concrete Structures - Part 1-1: General Rules for Buildings. DIN EN 1992-1:2011-1 (E).
展開 I-DEAS軟件J介紹!
仿真的廣泛應用:
5.非線性求解器（Model Solution-Non-Linear）：
支持幾何非線性、材料非線性、彈塑性及綜合非線性分析，利用Newton-Raphson方法求解非線性方程組。
6.變量化分析（Variationl Analysis）：
變量化分析仿真在設計初期介入設計過程，利用單一模型進行廣泛的設計研究，通過一次網格劃分和解算，生成手冊式結果，得到多種可對比的方案。
7.響應分析（Response Analysis）：
響應分析用來研究結構在靜態、瞬態、諧波和隨機等激勵下的受迫響應，模態可以來自結構分析或測試。
8.復合鋪合分析（Laminnate Composites）：
以復合鋪層材料結構進行高效的設計和評估.
9.注塑冷卻顧問（Part Advisor）：
注塑過程顧問系統。簡單實用，直接對STL格式進行計算，只用定義零件的材料和模具特性以及澆注口，就可以模擬澆注模具過程中塑料流動。可優化零件與模具的設計，以達到質量、成本和時間的最優平衡。
10.機械仿真（Mechanism Sim）：
分析機構在外力作用下的運動和受力。包含機構運動（Mechanism Design）的全部功能.
11.產品壽命預測（Durability）：
預測靜態載荷下產品的材料強度和疲勞安全的工具。
12.高級產品壽命預測（Advanced Durability）：
預測靜態或瞬態載荷下產品的產品壽命和疲勞破壞，包含產品壽命預測的全部功能。
13.電子系統冷卻仿真（Electronic System Cooling）：
電子系統的三維熱/流動分析。I-deas集成環境中的電子系統可靠性分析工具，仿真可直接在CAD模型上分析單元件、多芯片組、散熱片、PCB、功率模組或完整系統的熱/流動行為，含PCB Modeler和ECAD接口。
展開 
你想知道的CFD仿真分析流程在這里！ 附計算流體動力學分析：CFD軟件原理與應用下載
對于這些方程組，數學上已有相應的解法，如線性方程組可采用Guass消去法或Guass-Seidel迭代法求解，而對非線性方程組，可采用Newton-Raphson方法。
07
判斷解的收斂性
判斷計算殘差是否收斂
判斷進出口質量是否守恒
監測物理量變化
08
可視化輸出計算結果
流線圖
云圖
儲能集裝箱溫度云圖
潔凈室壓力云圖
車間速度云圖
矢量圖
車間速度矢量圖
空氣齡
潔凈室空氣齡圖
粒子沉積
下載地址：計算流體動力學分析：CFD軟件原理與應用
Abaqus應用之常見問題處理 （三）
2)ABAQUS/Standard使用Newton-Raphson方法來求解大型方程組，在每個增量步中都進行迭代運算，因此需要占用大量的硬盤和內存空間。其計算結果具有無條件穩定性，適合解決一般的線性和非線性問題。
3)ABAQUS/Explicit則采用中心差分法求解大型方程組，無需迭代計算，但需要進行大量的增量步計算，其計算時間取決于穩定極限值。其結果為有條件穩定，適用于高速動力學、復雜接觸、材料磨損和失效等問題的分析。
4)影響ABAQUS/Explicit穩定極限值的因素包括最小單元尺寸、材料特性、材料密度以及單元類型等。為了提高穩定極限值并加快分析速度，應選擇合適的網格密度，確保單元形狀規則，并在必要時應用質量縮放技術。
5)在ABAQUS/Explicit分析中，為避免數值振蕩，通常需要為模型定義阻尼，如體積粘性、材料阻尼、阻尼器單元和粘性壓力等。
6)對于ABAQUS/Explicit顯式求解器而言，模型中的位移邊界條件不應出現劇烈的變化，應該使用帶有平滑系數(Smoothing)的幅值曲線進行設置。
7) 靜態分析與動態分析的本質區別在于靜態分析不考慮慣性和沖擊效應，通常情況下，運動速度和加載速度對結果沒有影響；而動態分析則考慮慣性和沖擊效應，當載荷、位移、速度和加速度等在短時間內發生劇烈變化時（即加載速度過快），可能會導致顯著的局部變形。
8) ABAQUS/Explicit也可用于準靜態分析，關鍵在于合理設置加載速度、分析步時間和質量縮放系數等模型參數，以使結果盡可能接近靜態分析的結果。判斷是否為準靜態分析的一個重要標準是，分析過程中模型的動能一般不應超過內能的5%至10%。
展開 案例40-具有熱效應的形狀記憶合金（SMA）
分析和求解控制
使用大變形（NLGEOM，ON）和非對稱矩陣對全Newton-Raphson方法（NROPT，UNSYM）進行非線性靜態分析。整個解決方案分三個加載步驟完成，如圖所示：
比較了BEAM188和SOLID185模型的結果。
結果和討論
步驟1中載荷W拉伸的彈簧致動器如圖所示：
最大位移為43 mm，大于原始長度28 mm。
在步驟2中，在利用形狀記憶效應加熱后，彈簧致動器恢復到最大位移10 mm。變形處于馬氏體狀態以支撐重量，如下所示：
在步驟3中，在冷卻至250 K后，彈簧致動器拉伸回其原始長度：
以下是致動器的位移歷史：
位移歷史表明，BEAM188和SOLID185模型具有相似的結果。
然而，BEAM188模型的效率要高得多，需要大約一個小時才能完成。相比之下，SOLID185模型需要八個多小時才能完成。
建議
要執行涉及SMA的類似類型的分析，請考慮以下內容：
• 材料模型的剛度通常不對稱，但默認情況下在求解中使用對稱矩陣。如果在求解過程中出現收斂困難，請指定非對稱求解器選項（NROPT、UNSYM）。
• SMA相變包括相變階段和飽和相變。因為在轉變階段更難實現收斂，所以將轉變載荷步分解為更小的載荷步，并使用更小的時間步長。
• 材料模型（TB、SMA）的超彈性選項（TBOPT=SUPE）僅支持三維、平面應變和軸對稱應力狀態。材料模型的記憶選項（TBOPT=MEFF）支持大多數應力狀態，包括梁、殼、平面應變、軸對稱和三維應力狀態。兩種材質選項都不支持LINK180單元。
參考文獻
Paremer, A., Fumer, S., Rice, D. P. (1992).
展開 ——科研狗的義務和擔當
求解方法常用的是Newton——Raphson 方法。怎么求K 怎么寫F，怎么編程，就是需要的。不同的人的角度，寫公式，編程的角度不同，導致了人為的困難。 我先建議看國外的寫得好的程序。
先在這拋一塊磚，什么時候有時間了，我去別的地方搬一塊玉給大家看看
有限元理論基礎及Abaqus內部實現方式研究系列4：非線性問題的求解 ￥1
2.1.2.3 Modified Newton-Raphson法
如果非線性曲線比較平緩，那么在New-Raphson迭代過程中迭代步的切線剛度矩陣可以近似為上幾個或者第一個迭代步的，即。這樣可以只在第一個迭代步計算一次K，提高每個迭代步的計算速度，但由于一般情況下
的變化較大，所以會降低收斂速度。
2.1.2.4 BFGS法
其實就是Abaqus非線性Step設置中的Quasi-Newton方法，在Abaqus做非線性可以在Newton方法和Quasi-Newton二選一。這里不做更多說明。
2.1.3 增量迭代法
一般情況增量法可以保證求解過程的收斂性但收斂速度較慢，而Newton-Raphson法收斂速度較慢但收斂性沒有保證，所以混合法結合了增量法和Newton-Raphson法來求解非線性問題。混合法首先將外載荷分成若干個增量步，在每個增量步內采用Newton-Raphson法迭代求解，在增量步內求解完成后繼續求解下一個增量步，最后將所有增量步累加起來即得到結果。
2.1.4 收斂判據
常見的收斂判據有分為兩種，失衡力準則、位移準則
2.1.4.1 失衡力準則
2.1.4.2 位移準則
2.2 Abaqus的非線性問題求解
Abaqus的非線性中上面說的三種迭代都涉及，其中General Static分析步采用Newton或者Quasi-Newton方法，而Static, Riks分析步采用Riks弧長法。我們現在只聚焦到General Static的Newton方法。
展開 轉，穩態傳熱分析
Command: NCNV
GUI: Main Menu>Solution>-Load Step Opts-Nolinear>Criteria to Stop
· 線性搜索：設置本選項可使ANSYS用Newton-Raphson方法進行線性搜索。
Command: LNSRCH
GUI: Main Menu>Solution>-Load Step Opts-Nolinear>Line Search
· 預測矯正：本選項可激活每一子步第一次迭代對自由度求解的預測矯正。
Command: PRED
GUI: Main Menu>Solution>-Load Step Opts-Nolinear>Predictor
c. 輸出控制
· 控制打印輸出：本選項可將任何結果數據輸出到*.out 文件中。
Command: OUTPR
GUI: Main Menu>Solution>-Load Step Opts-Output Ctrls>Solu Printout
· 控制結果文件：控制*.rth的內容。
Command: OUTRES
GUI: Main Menu>Solution>-Load Step Opts-Output Ctrls>DB/Results File
④、確定分析選項
a. Newton-Raphson選項（僅對非線性分析有用）
Command: NROPT
GUI: Main Menu>Solution>Analysis Options
b.
展開 【接觸分析】詳述ABAQUS接觸分析（1/2）
4.5 接觸算法
ABAQUS/Standard中的接觸算法是圍繞Newton-Raphson方法建立的，在每個增量步開始時檢查所有接觸相互作用的狀態，以建立從屬節點是開放還是閉合的。圖中，
p表示從屬節點上的接觸壓力，
h表示從屬節點侵入主控表面的距離。如果一個節點是閉合的，ABAQUS/Standard需要確定它是處于滑動還是粘結。ABAQUS/Standard對每個閉合節點施加一個約束，而對那些改變接觸狀態從閉合到開放的任何節點接觸約束。然后，ABAQUS/Standard進行迭代，并利用計算的修正值來更新模型的構形。
圖為接觸算法流程圖
ABAQUS/Standard首先檢驗在從屬節點上接觸條件的變化，如果在當前迭代步中檢測到任何接觸變化，ABAQUS/Standard將標識其為 嚴重不連續迭代（severe discontinuity iteration），并不再進行力或力矩的平衡檢驗。在第1次迭代后，ABAQUS/Standard修正接觸約束以反映接觸狀態的變化，并試圖進行第2次迭代，重復此過程，直到完成迭代并且不改變接觸狀態。
對于每個完成的增量步，在信息文件和狀態文件中的總結將顯示出有多少次迭代是嚴重不連續迭代和多少次是平衡迭代。對于每個增量步，其總的迭代數目是這兩者之和。如果嚴重不連續迭代的數目很多，并只有很少的平衡迭代，那么ABAQUS/Standard難以確定合適的接觸狀態。對于任何需要嚴重不連續迭代超過 12次的增量步，ABAQUS/Standard會默認放棄，而應用更小的增量尺度再次進行增量步計算。如果沒有嚴重不連續迭代，從一個增量步到下一個增量步的接觸狀態不會發生改變。
展開 
案例52-粘土層上路堤的順序施工
• 使用非對稱Newton-Raphson解方法（NROPT，UNSYM）和耦合孔隙壓力熱機械單元（本問題中為CPT212）。
• 忽略可能在地質靜態求解開始時出現的接觸單元狀態警告。（由于接觸單元被殺死，然后被重新激活，預計狀態會突然改變。）
在 COMSOL 中如何求解積分上下限
由于我們使用 Newton-Raphson 方法來求解這個問題，因此不應該從函數斜率為零的點開始。經過求解，我們發現 u_b = 2.621。
現在，我們讓事情變得復雜一點，這次需要求解的是
和
這兩個極限：
由于我們有兩個未知數，顯然需要在這里再加一個方程，所以讓我們再假設
。
添加了一個額外的方程來指定區間上限和下限之間的差值。
求解如上所示的模型，將得到 u_a = 1.932 和 u_b = 2.932 的值。實際上也可以用單個全局方程來求解這個問題，方法是將 6-integrate(u^2, u, u_b-1, u_b) 寫成求解 u_b 的方程，但有趣的是，我們可以同時求解多個方程。
一個傳熱示例
接下來，讓我們將上文中介紹的技術用于實踐，來確定換熱器的運行條件。以 COMSOL 案例庫中的一個示例來說明，即通過一個浸沒在池塘中的管道網絡循環對水進行地熱加熱。
通過水下管道網絡泵送的水被加熱。
在這個示例中，管道流模塊用于模擬 5°C(278.15 K)的水被泵入管網，并被池塘中相對較溫暖的水加熱。池塘中的水溫隨深度在 10°C 至 20°C 之間變化，輸出端的計算溫度為 11.1°C(284.25 K)。如果水的質量流量指定為 4 kg/s，那么吸收的總熱量為：
其中，
是質量流量，
是與溫度有關的比熱。
實際上，這里的管道網絡是一個閉環系統，但我們根本沒有對管道出口和入口之間的系統部分進行建模，這里包含一個隱含的假設，就是當水從出口泵回入口時，它被精確冷卻到 5°C。
因此，與其假設進入管道的水的溫度是恒定溫度，不如考慮這個閉環系統連接到另一個換熱器，通過它來去除指定的熱量。
展開 『轉貼』沖壓件坯料的優化設計方法研究
3.3 對塑性變形功求極值
為了求整體塑性功的極小值，并且考慮壓邊力、拉延筋等邊界條件，令：
(9)
求的最小值，且令
(10)
用Newton-Raphson 方法解非線性方程組（1），即得：
(11)
(12)
其中是減速因子,取值為0到1。
3.4 一步模擬的實例：盒形件拉深 [21]
圖1模具尺寸圖示
Fig.1 Tooling geometry for the deep drawing of square cup.
盒形件及模具尺寸如圖1所示。材料的性能指數及工藝參數為：應力應變關系曲線為 （MPa）；厚向異性參數；板料厚度（mm）；摩擦系數。圖2表示的是用三角形單元離散化后的工件。
利用一步模擬算法，在微機（Windows2000, CPU Intel P4/2GHz）上計算僅需要2min，計算效率大大高于增量法有限元。可以快速獲取優化板料形狀（圖3所示），并且可以預測最終工件的厚向應力應變分布。由圖 4可知，利用該優化的坯料，可以獲得較好的拉深成形性能。
展開 ABAQUS子程序UMAT里彈塑本構的實現
另外，每一次求解都需要對求倒數，如果求解方程組，就是對剛度矩陣求逆，這種方法在求解中控制兩次求解之差，當其值很小時，就認為接近真實值了，迭代結束
圖4-3
2：Newton-Raphson方法
Newton-Raphson方法的算法與常剛度法不同[12]，如果得近似表達式是不成立的，存在著殘余值，即，此式也可以作為近似值與真實值的差值量度，實際上在具體計算時，也可以控制其值，當極小時，就認為接近真實值了，當第次迭代的值
是真實解，則可以按照Taylor級數展開得到
圖4-4
3：切線剛度法
在復雜非線性問題求解中，剛度與的大小是有一定關系的，在用增量法來求解這種問題時，就等于結構任一點處力與位移的曲線的局部梯度，稱為切線剛度[13]，剛度矩陣的倒數很難用自變量顯示表達，通過增量方式求解，在每一步荷載增量范圍內把問題線性化，求解方法與Newton-Raphson方法相同
總結以上可以得到：
以上幾種算法中，通過比較，不難發現，直接迭代法采用了固定的剛度，適合解決非線性程度不高的本構關系，而切線剛度法采用了變化的剛度，在每一步上都做了實時的修正，對非線性程度較高本構關系任然有效，在效率和迭代精度方面，切線剛度法采用的修正更符合非線性材料的應力應變關系，具有較大的優勢，這也是本文采用切線剛度法計算的原因，當然，非線性有限元算法還有很多，切線剛度法也不見得就是最好的能解決所有問題的算法，但是它是在程序開發難度不高和精度方面較高的條件下相對來說最好的
本文采用的本構關系是同性硬化彈塑性模型[14]，采用Mises屈服準則，下面將根據J2理論[15]，分別推導常剛度法和切線剛度法計算該問題的的算法公式
展開 