本發(fā)明涉及數(shù)值計(jì)算領(lǐng)域，具體而言，涉及一種基于傳輸線迭代法的非線性靜磁場模型的有限元求解方法，該方法主要針對(duì)2D軸對(duì)稱非線性靜態(tài)電磁場進(jìn)行求解。

背景技術(shù)：

有限元法是工業(yè)設(shè)計(jì)中最常用的數(shù)值計(jì)算方法，被諸多商用仿真軟件采用，應(yīng)用廣泛。然而，隨著求解模型的日益復(fù)雜化以及分網(wǎng)單元數(shù)目的不斷增多，以傳統(tǒng)的牛頓迭代法為核心的非線性有限元求解方法面臨著求解耗時(shí)嚴(yán)重的問題，這直接關(guān)系到產(chǎn)品研發(fā)的速度和效率。

有限元問題的求解的核心在于求解線性方程組，而對(duì)于非線性問題來說，傳統(tǒng)的牛頓迭代法每一步都要利用新的迭代結(jié)果重新生成有限元模型的全局矩陣，隨著模型分網(wǎng)的不斷增大，全局矩陣的維度不斷變大，每一步矩陣的LU分解等消耗的時(shí)間會(huì)相應(yīng)的增大，總體的求解時(shí)間可能隨著分網(wǎng)的變密而成幾何式增大。

因此，需要研究一種新的迭代方法，以解決牛頓迭代法求解有限元非線性問題時(shí)帶來的求解時(shí)間長，效率低的問題。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明提供了一種基于傳輸線迭代法的2D軸對(duì)稱非線性靜磁場模型的有限元求解方法，用以解決牛頓迭代法求解有限元非線性問題時(shí)帶來求解時(shí)間長，效率低的問題。

為了達(dá)到上述目的，本發(fā)明提供了一種基于傳輸線迭代法的2D軸對(duì)稱非線性靜磁場模型的有限元求解方法，其包括以下步驟：

S1：確定待求解的變量以及求解域，待求解的變量為一2D軸對(duì)稱非線性靜磁場的磁勢(shì)，2D軸對(duì)稱非線性靜磁場由通電線圈中的電流產(chǎn)生，通電線圈周圍的各元件均為鐵磁材料，求解域?yàn)?D軸對(duì)稱非線性靜磁場所在的區(qū)域；

S2：建立一平面x-y坐標(biāo)系，以2D軸對(duì)稱非線性靜磁場的對(duì)稱軸為y軸，在y軸上選定其中一點(diǎn)為原點(diǎn)，并設(shè)定經(jīng)過原點(diǎn)并與y軸垂直的直線為x軸，即x-y平面為2D軸對(duì)稱非線性靜磁場所在區(qū)域一過對(duì)稱軸的截面所在的平面；

S3：列出2D軸對(duì)稱非線性靜磁場中的控制方程和邊界條件式并組成一微分方程組，其控制方程為：

其中，J為電流密度變量，μ為三角單元的磁導(dǎo)率，A為磁勢(shì)，

邊界條件式為：

Γ₁：A＝0，

Γ₁表示磁勢(shì)A在邊界Γ₁上的分布，Γ₂表示磁勢(shì)A沿邊界的外法線方向的變化率，

S4：使用三角單元對(duì)求解域進(jìn)行分網(wǎng)，得到包含多個(gè)三角單元的有限元網(wǎng)絡(luò)，該有限元網(wǎng)絡(luò)中的三角單元總個(gè)數(shù)為N，節(jié)點(diǎn)總個(gè)數(shù)為M，并分別對(duì)三角單元和節(jié)點(diǎn)進(jìn)行1～N和1～M的編號(hào)，其中1000≤N≤3000；

S5：根據(jù)微分方程組的泛函形式，推導(dǎo)出每一個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]，其中，每一[Y^e]均為3×3的矩陣，每一[J^e]均為1×3的矩陣：

[J^e]＝[J_l J_m J_n]，

l、m、n分別為推導(dǎo)每一三角單元的[Y^e]和[J^e]時(shí)，三角單元的三個(gè)頂點(diǎn)的編號(hào)，

r和s分別為三角單元的三個(gè)頂點(diǎn)編號(hào)1、m和n中的其中兩個(gè)頂點(diǎn)編號(hào)，

x₁、x_m和x_n分別為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n在平面坐標(biāo)系中的橫坐標(biāo)，y₁、y_m和y_n分別為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n在平面坐標(biāo)系中的縱坐標(biāo)，Δ為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n組成的三角單元的面積；

S6：根據(jù)得到的每一個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]，對(duì)N個(gè)三角單元進(jìn)行有限元裝配，得到全局矩陣Y和J，其中Y為M×M矩陣，J為M×1矩陣；

S7：求解非線性方程組YA＝J，得到2D軸對(duì)稱非線性靜磁場中每個(gè)節(jié)點(diǎn)的磁勢(shì)A，其中A為M×1的節(jié)點(diǎn)磁勢(shì)矩陣，A＝[A₁ A₂ … A_M]^T；

S8：根據(jù)步驟S8中計(jì)算得到的節(jié)點(diǎn)磁勢(shì)矩陣A，按照以下各式計(jì)算每一個(gè)三角單元的磁感應(yīng)強(qiáng)度B，其中，

S9：根據(jù)鐵磁材料的B-H曲線以及步驟S8中計(jì)算得到的每一個(gè)三角單元的磁感應(yīng)強(qiáng)度B，并計(jì)算出每一個(gè)三角單元的磁導(dǎo)率μ；

S10：以步驟S4中的分網(wǎng)結(jié)果為基礎(chǔ)，對(duì)求解域進(jìn)行精細(xì)的三角分網(wǎng)，得到三角單元總個(gè)數(shù)為N'、節(jié)點(diǎn)總個(gè)數(shù)為M'的有限元網(wǎng)絡(luò)，并分別對(duì)三角單元和節(jié)點(diǎn)進(jìn)行1～N'和1～M'的編號(hào)；

S11：按照步驟S5中的方法，對(duì)步驟S10中得到的有限元網(wǎng)絡(luò)再次計(jì)算每個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]；

S12：有限元網(wǎng)絡(luò)轉(zhuǎn)化為電路模型，將步驟S11中得到的單元矩陣[Y^e]視為電路的導(dǎo)納矩陣，激勵(lì)源單元矩陣[J^e]視為每個(gè)節(jié)點(diǎn)與地之間的電流源矩陣，對(duì)有限元網(wǎng)絡(luò)中的每一個(gè)三角單元均建立一個(gè)等效電路網(wǎng)絡(luò)，建立等效電路網(wǎng)絡(luò)的方法如下：

將單元矩陣[Y^e]對(duì)角線上的元素視為自導(dǎo)，非對(duì)角線上的元素視為互導(dǎo)，

對(duì)于非對(duì)角線上的元素，若Y_rs＞0，則在三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)中的節(jié)點(diǎn)r和節(jié)點(diǎn)s之間設(shè)置一受控電流源，該受控電流源中的電流大小為U_rsY_rs,方向?yàn)閺墓?jié)點(diǎn)r流向節(jié)點(diǎn)s，其中U_rs為節(jié)點(diǎn)r與節(jié)點(diǎn)s之間的磁勢(shì)差，

對(duì)于非對(duì)角線上的元素，若Y_rs＜0，則在三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)中的節(jié)點(diǎn)r和節(jié)點(diǎn)s之間設(shè)置一純電阻，該純電阻的導(dǎo)納為|Y_rs|，

若有限元單元矩陣[Y^e]的第r行的所有元素之和不等于0，當(dāng)?shù)趓行所有元素之和大于零時(shí)，在節(jié)點(diǎn)r與地之間設(shè)置一純電阻，該純電阻的導(dǎo)納為Y_rl+Y_rm+Y_rn，當(dāng)?shù)趓行所有元素之和小于零時(shí)，則在節(jié)點(diǎn)r與地之間設(shè)置一受控電流源，該受控電流源中的電流大小為U_r0·|Y_rl+Y_rm+Y_rn|,方向?yàn)閺墓?jié)點(diǎn)r流向地，其中U_r0為節(jié)點(diǎn)r與地之間的磁勢(shì)差，

在每一節(jié)點(diǎn)與地之間均設(shè)置一電流源，節(jié)點(diǎn)l、節(jié)點(diǎn)m、節(jié)點(diǎn)n與地之間的電流源中的電流大小分別為J_l、J_m、J_n，電流方向?yàn)橛傻亓飨蚬?jié)點(diǎn)；

S13：組裝電路，將步驟S12中建立的每一個(gè)三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)通過節(jié)點(diǎn)進(jìn)行連接，組裝成一個(gè)完整的非線性電路網(wǎng)絡(luò)，該非線性電路網(wǎng)絡(luò)等效為包含一線性網(wǎng)絡(luò)與多個(gè)非線性待求元件的電路；

S14：對(duì)于步驟S13中得到的非線性電路網(wǎng)絡(luò)，為了使用傳輸線迭代方法進(jìn)行求解，需要在非線性元件與線性網(wǎng)絡(luò)之間添加一段傳輸線，由于傳輸線對(duì)信號(hào)傳輸?shù)难訒r(shí)作用，電路的非線性求解過程包括入射階段和反射階段，

入射階段，非線性電路元件的電壓信號(hào)向線性網(wǎng)絡(luò)進(jìn)行入射，等效為傳輸線導(dǎo)納與虛擬電流源的并聯(lián)，

反射階段，電壓信號(hào)由線性網(wǎng)絡(luò)傳向非線性元件，對(duì)非線性元件進(jìn)行求解，如此不斷迭代入射階段和反射階段，直到電路達(dá)到穩(wěn)態(tài)，

(一)在線性部分與非線性元件之間添加一段傳輸線，傳輸線的導(dǎo)納的計(jì)算方法如下：

(1)確定每一個(gè)三角單元的磁導(dǎo)率μ的估計(jì)值，檢查經(jīng)過步驟S10分網(wǎng)后得到的三角單元的重心對(duì)應(yīng)的第一次分網(wǎng)的三角單元，并將對(duì)應(yīng)的第一次分網(wǎng)的三角單元的磁導(dǎo)率設(shè)為三角單元的磁導(dǎo)率，

(2)非線性元件的導(dǎo)納是一個(gè)關(guān)于磁導(dǎo)率μ的變量，將上一步得到的μ值代入到非線性元件表達(dá)式，得到的結(jié)果作為對(duì)應(yīng)的傳輸線的導(dǎo)納值，

(二)設(shè)迭代開始時(shí)每一節(jié)點(diǎn)的電壓均為0，當(dāng)?shù)趎個(gè)節(jié)點(diǎn)電壓信號(hào)以V_in反射到線性網(wǎng)絡(luò)時(shí)，每一非線性待求元件等效為一導(dǎo)納和一電流源的并聯(lián)電路，其中，導(dǎo)納為對(duì)應(yīng)的傳輸線導(dǎo)納Y_n，電流源中的電流值為2V_inY_n，對(duì)該等效電路進(jìn)行求解，得到每一節(jié)點(diǎn)的電壓值V_in，

(三)根據(jù)各個(gè)節(jié)點(diǎn)的電壓值，利用非線性元件與電壓之間的關(guān)系式，計(jì)算并更新非線性元件的導(dǎo)納值，

(四)計(jì)算各個(gè)節(jié)點(diǎn)向非線性元件入射的電壓值V_rn，節(jié)點(diǎn)n處的V_rn＝V_n-V_in，

(五)入射過程，每一非線性待求元件等效為一導(dǎo)納與一電流源的并聯(lián)電路，其中，導(dǎo)納為對(duì)應(yīng)的傳輸線導(dǎo)納Y_n，電流源中的電流值為2V_rnY_n，得到每一非線性元件兩端的電壓

(六)計(jì)算各個(gè)節(jié)點(diǎn)向線性網(wǎng)絡(luò)入射的電壓值V_in，節(jié)點(diǎn)n處的

(七)重復(fù)步驟(二)～(六)，直至相鄰兩次迭代中，步驟(二)所求的電壓值V_n達(dá)到預(yù)設(shè)的收斂誤差，此時(shí)計(jì)算得到的各節(jié)點(diǎn)的電壓值V_n即為所求電壓值，

S15：根據(jù)每一節(jié)點(diǎn)的電壓值繪制2D軸對(duì)稱非線性靜磁場中的磁勢(shì)云圖。

在本發(fā)明的一實(shí)施例中，步驟S14中(二)對(duì)等效電路的求解為使用節(jié)點(diǎn)電壓法進(jìn)行求解，其步驟為：

(一)計(jì)算得到矩陣YV＝I，其中Y為電路導(dǎo)納矩陣，由于迭代過程中，導(dǎo)納矩陣Y保持不變，僅需計(jì)算一次，V為待求節(jié)點(diǎn)電壓，I為節(jié)點(diǎn)電流；

(二)迭代第一步對(duì)矩陣Y進(jìn)行LU分解，即Y＝LU，其中L為單位下三角矩陣，U為上三角矩陣，之后的每一次迭代，無需計(jì)算此步，直接計(jì)算步驟(三)；

(三)使用公式V＝U^-1(L^-1I)求解節(jié)點(diǎn)電壓V。

在本發(fā)明的一實(shí)施例中，步驟S14中(五)，入射過程中，每一非線性元件兩端電壓的求解是獨(dú)立的，在此使用多核并行技術(shù)同時(shí)對(duì)多個(gè)非線性元件兩端電壓進(jìn)行求解。

本發(fā)明提供的基于傳輸線迭代法的2D軸對(duì)稱非線性靜磁場模型的有限元求解方法解決了牛頓迭代法求解有限元非線性問題時(shí)帶來求解時(shí)間長，效率低的問題。本發(fā)明中的每一步驟均無需再次計(jì)算全局矩陣，只需要進(jìn)行一次全局矩陣的LU分解，然后即可重復(fù)使用，從而節(jié)省計(jì)算時(shí)間；同時(shí)，該方法非常適合使用并行算法進(jìn)行加速，能夠進(jìn)一步的加快有限元問題的求解。相對(duì)于傳統(tǒng)的牛頓迭代法，本發(fā)明在求解時(shí)間上有著很大的優(yōu)勢(shì)，有著廣闊的應(yīng)用前景。

附圖說明

為了更清楚地說明本發(fā)明實(shí)施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對(duì)實(shí)施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實(shí)施例，對(duì)于本領(lǐng)域普通技術(shù)人員來講，在不付出創(chuàng)造性勞動(dòng)的前提下，還可以根據(jù)這些附圖獲得其他的附圖。

圖1為產(chǎn)生一2D軸對(duì)稱非線性靜磁場的接觸器機(jī)械結(jié)構(gòu)示意圖；

圖2為與圖1對(duì)應(yīng)的靜磁場的有限元模型區(qū)域示意圖；

圖3a為對(duì)求解域進(jìn)行首次分網(wǎng)(粗糙分網(wǎng))的示意圖；

圖3b為對(duì)求解域進(jìn)行再次分網(wǎng)(精細(xì)分網(wǎng))的示意圖；

圖4為求解域中的三角單元的示意圖；

圖5為等效電路網(wǎng)絡(luò)的示意圖；

圖6為對(duì)三角單元進(jìn)行組裝的示意圖；

圖7為傳輸線迭代等效示意圖；

圖8為反射過程的等效示意圖；

圖9為入射過程的等效示意圖；

圖10為靜磁場中的磁勢(shì)云圖；

圖11為傳統(tǒng)牛頓迭代法與傳輸線迭代法求解時(shí)間對(duì)比；

圖12為牛頓迭代法與傳輸線迭代法不同分網(wǎng)大小計(jì)算時(shí)間對(duì)比；

圖13為牛頓迭代法與傳輸線迭代法單步計(jì)算時(shí)間對(duì)比。

具體實(shí)施方式

下面將結(jié)合本發(fā)明實(shí)施例中的附圖，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚、完整地描述，顯然，所描述的實(shí)施例僅僅是本發(fā)明一部分實(shí)施例，而不是全部的實(shí)施例?；诒景l(fā)明中的實(shí)施例，本領(lǐng)域普通技術(shù)人員在沒有付出創(chuàng)造性勞動(dòng)前提下所獲得的所有其他實(shí)施例，都屬于本發(fā)明保護(hù)的范圍。

本發(fā)明提供的基于傳輸線迭代法的2D軸對(duì)稱非線性靜磁場模型的有限元求解方法包括以下步驟：

S1：確定待求解的變量以及求解域，待求解的變量為一2D軸對(duì)稱非線性靜磁場的磁勢(shì)，2D軸對(duì)稱非線性靜磁場由通電線圈中的電流產(chǎn)生，通電線圈周圍的各元件均為鐵磁材料，求解域?yàn)?D軸對(duì)稱非線性靜磁場所在的區(qū)域；

圖1為產(chǎn)生一2D軸對(duì)稱非線性靜磁場的接觸器機(jī)械結(jié)構(gòu)示意圖，如圖所示，圖1中為一接觸器，其中的線圈中通有電流，線圈周圍的鐵芯、推桿、銜鐵等的材質(zhì)均為鐵磁材料，該接觸器周圍的靜磁場以推桿為中心而軸對(duì)稱，因此，可以先以靜磁場的其中任意一個(gè)截面(如圖中的方框部分)為研究對(duì)象。

圖2為與圖1對(duì)應(yīng)的靜磁場的有限元模型區(qū)域示意圖，該有限元模型區(qū)域?yàn)閳D1中的俯視區(qū)域的右側(cè)部分，該區(qū)域也就是本發(fā)明針對(duì)的求解域。

S2：建立一平面x-y坐標(biāo)系，以2D軸對(duì)稱非線性靜磁場的對(duì)稱軸為y軸，在y軸上選定其中一點(diǎn)為原點(diǎn)，并設(shè)定經(jīng)過原點(diǎn)并與y軸垂直的直線為x軸，即x-y平面為2D軸對(duì)稱非線性靜磁場所在區(qū)域一過對(duì)稱軸的截面所在的平面；

S3：列出2D軸對(duì)稱非線性靜磁場中的控制方程和邊界條件式并組成一微分方程組，其控制方程為：

其中，J為電流密度變量，μ為三角單元的磁導(dǎo)率，A為磁勢(shì)，

邊界條件式為：

Γ₁：A＝0，

Γ₁表示磁勢(shì)A在邊界Γ₁上的分布，Γ₂表示磁勢(shì)A沿邊界的外法線方向的變化率，

S4：使用三角單元對(duì)求解域進(jìn)行分網(wǎng)，得到包含多個(gè)三角單元的有限元網(wǎng)絡(luò)，該有限元網(wǎng)絡(luò)中的三角單元總個(gè)數(shù)為N，節(jié)點(diǎn)總個(gè)數(shù)為M，并分別對(duì)三角單元和節(jié)點(diǎn)進(jìn)行1～N和1～M的編號(hào)，其中1000≤N≤3000；

圖3為對(duì)求解域進(jìn)行首次分網(wǎng)(粗糙分網(wǎng))的示意圖。

S5：根據(jù)微分方程組的泛函形式，推導(dǎo)出每一個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]，其中，每一[Y^e]均為3×3的矩陣，每一[J^e]均為1×3的矩陣：

[J^e]＝[J_l J_m J_n]，

圖4為求解域中的三角單元的示意圖，l、m、n分別為推導(dǎo)每一三角單元的[Y^e]和[J^e]時(shí)，三角單元的三個(gè)頂點(diǎn)的編號(hào)，

r和s分別為三角單元的三個(gè)頂點(diǎn)編號(hào)1、m和n中的其中兩個(gè)頂點(diǎn)編號(hào)，

x₁、x_m和x_n分別為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n在平面坐標(biāo)系中的橫坐標(biāo)，y¹、y_m和y_n分別為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n在平面坐標(biāo)系中的縱坐標(biāo)，Δ為節(jié)點(diǎn)l、節(jié)點(diǎn)m和節(jié)點(diǎn)n組成的三角單元的面積；

S6：根據(jù)得到的每一個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]，對(duì)N個(gè)三角單元進(jìn)行有限元裝配，得到全局矩陣Y和J，其中Y為M×M矩陣，J為M×1矩陣；

S7：求解非線性方程組YA＝J，得到2D軸對(duì)稱非線性靜磁場中每個(gè)節(jié)點(diǎn)的磁勢(shì)A，其中A為M×1的節(jié)點(diǎn)磁勢(shì)矩陣，A＝[A₁ A₂…A_M]^T；

S8：根據(jù)步驟S8中計(jì)算得到的節(jié)點(diǎn)磁勢(shì)矩陣A，按照以下各式計(jì)算每一個(gè)三角單元的磁感應(yīng)強(qiáng)度B，其中，

S9：根據(jù)鐵磁材料的B-H曲線以及步驟S8中計(jì)算得到的每一個(gè)三角單元的磁感應(yīng)強(qiáng)度B，并計(jì)算出每一個(gè)三角單元的磁導(dǎo)率μ；

S10：以步驟S4中的分網(wǎng)結(jié)果為基礎(chǔ)，對(duì)求解域進(jìn)行精細(xì)的三角分網(wǎng)，圖3b為對(duì)求解域進(jìn)行再次分網(wǎng)(精細(xì)分網(wǎng))的示意圖，得到三角單元總個(gè)數(shù)為N'、節(jié)點(diǎn)總個(gè)數(shù)為M'的有限元網(wǎng)絡(luò)，并分別對(duì)三角單元和節(jié)點(diǎn)進(jìn)行1～N'和1～M'的編號(hào)；

S11：按照步驟S5中的方法，對(duì)步驟S10中得到的有限元網(wǎng)絡(luò)再次計(jì)算每個(gè)三角單元的單元矩陣[Y^e]和激勵(lì)源單元矩陣[J^e]；

S12：有限元網(wǎng)絡(luò)轉(zhuǎn)化為電路模型，將步驟S11中得到的單元矩陣[Y^e]視為電路的導(dǎo)納矩陣，激勵(lì)源單元矩陣[J^e]視為每個(gè)節(jié)點(diǎn)與地之間的電流源矩陣，對(duì)有限元網(wǎng)絡(luò)中的每一個(gè)三角單元均建立一個(gè)等效電路網(wǎng)絡(luò)，圖5為等效電路網(wǎng)絡(luò)的示意圖，如圖所示，建立等效電路網(wǎng)絡(luò)的方法如下：

將單元矩陣[Y^e]對(duì)角線上的元素視為自導(dǎo)，非對(duì)角線上的元素視為互導(dǎo)，

對(duì)于非對(duì)角線上的元素，若Y_rs＞0，則在三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)中的節(jié)點(diǎn)r和節(jié)點(diǎn)s之間設(shè)置一受控電流源，該受控電流源中的電流大小為U_rsY_rs,方向?yàn)閺墓?jié)點(diǎn)r流向節(jié)點(diǎn)s，其中U_rs為節(jié)點(diǎn)r與節(jié)點(diǎn)s之間的磁勢(shì)差，

對(duì)于非對(duì)角線上的元素，若Y_rs＜0，則在三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)中的節(jié)點(diǎn)r和節(jié)點(diǎn)s之間設(shè)置一純電阻，該純電阻的導(dǎo)納為|Y_rs|，

若有限元單元矩陣[Y^e]的第r行的所有元素之和不等于0，當(dāng)?shù)趓行所有元素之和大于零時(shí)，在節(jié)點(diǎn)r與地之間設(shè)置一純電阻，該純電阻的導(dǎo)納為Y_rl+Y_rm+Y_rn，當(dāng)?shù)趓行所有元素之和小于零時(shí)，則在節(jié)點(diǎn)r與地之間設(shè)置一受控電流源，該受控電流源中的電流大小為U_r0·|Y_rl+Y_rm+Y_rn|,方向?yàn)閺墓?jié)點(diǎn)r流向地，其中U_r0為節(jié)點(diǎn)r與地之間的磁勢(shì)差，

在每一節(jié)點(diǎn)與地之間均設(shè)置一電流源，節(jié)點(diǎn)l、節(jié)點(diǎn)m、節(jié)點(diǎn)n與地之間的電流源中的電流大小分別為J_l、J_m、J_n，電流方向?yàn)橛傻亓飨蚬?jié)點(diǎn)；

S13：組裝電路，圖6為對(duì)三角單元進(jìn)行組裝的示意圖，如圖所示，將步驟S12中建立的每一個(gè)三角單元對(duì)應(yīng)的等效電路網(wǎng)絡(luò)通過節(jié)點(diǎn)進(jìn)行連接，組裝成一個(gè)完整的非線性電路網(wǎng)絡(luò)，該非線性電路網(wǎng)絡(luò)等效為包含一線性網(wǎng)絡(luò)與多個(gè)非線性待求元件的電路；

S14：對(duì)于步驟S13中得到的非線性電路網(wǎng)絡(luò)，為了使用傳輸線迭代方法進(jìn)行求解，需要在非線性元件與線性網(wǎng)絡(luò)之間添加一段傳輸線，由于傳輸線對(duì)信號(hào)傳輸?shù)难訒r(shí)作用，電路的非線性求解過程包括入射階段和反射階段，圖7為傳輸線迭代等效示意圖，

入射階段，非線性電路元件的電壓信號(hào)向線性網(wǎng)絡(luò)進(jìn)行入射，等效為傳輸線導(dǎo)納與虛擬電流源的并聯(lián)，

反射階段，電壓信號(hào)由線性網(wǎng)絡(luò)傳向非線性元件，對(duì)非線性元件進(jìn)行求解，如此不斷迭代入射階段和反射階段，直到電路達(dá)到穩(wěn)態(tài)，

(一)在線性部分與非線性元件之間添加一段傳輸線，傳輸線的導(dǎo)納的計(jì)算方法如下：

(1)確定每一個(gè)三角單元的磁導(dǎo)率μ的估計(jì)值，檢查經(jīng)過步驟S10分網(wǎng)后得到的三角單元的重心對(duì)應(yīng)的第一次分網(wǎng)的三角單元，并將對(duì)應(yīng)的第一次分網(wǎng)的三角單元的磁導(dǎo)率設(shè)為三角單元的磁導(dǎo)率，

(2)非線性元件的導(dǎo)納是一個(gè)關(guān)于磁導(dǎo)率μ的變量，將上一步得到的μ值代入到非線性元件表達(dá)式，得到的結(jié)果作為對(duì)應(yīng)的傳輸線的導(dǎo)納值，

(二)設(shè)迭代開始時(shí)每一節(jié)點(diǎn)的電壓均為0，當(dāng)?shù)趎個(gè)節(jié)點(diǎn)電壓信號(hào)以V_in反射到線性網(wǎng)絡(luò)時(shí)，每一非線性待求元件等效為一導(dǎo)納和一電流源的并聯(lián)電路，其中，導(dǎo)納為對(duì)應(yīng)的傳輸線導(dǎo)納Y_n，電流源中的電流值為2V_inY_n，對(duì)該等效電路進(jìn)行求解，得到每一節(jié)點(diǎn)的電壓值V_in，圖8為反射過程的等效示意圖，

(三)根據(jù)各個(gè)節(jié)點(diǎn)的電壓值，利用非線性元件與電壓之間的關(guān)系式，計(jì)算并更新非線性元件的導(dǎo)納值，

(四)計(jì)算各個(gè)節(jié)點(diǎn)向非線性元件入射的電壓值V_rn，節(jié)點(diǎn)n處的V_rn＝V_n-V_in，

(五)入射過程，每一非線性待求元件等效為一導(dǎo)納與一電流源的并聯(lián)電路，其中，導(dǎo)納為對(duì)應(yīng)的傳輸線導(dǎo)納Y_n，電流源中的電流值為2V_rnY_n，得到每一非線性元件兩端的電壓圖9為入射過程的等效示意圖，

(六)計(jì)算各個(gè)節(jié)點(diǎn)向線性網(wǎng)絡(luò)入射的電壓值V_in，節(jié)點(diǎn)n處的

(七)重復(fù)步驟(二)～(六)，直至相鄰兩次迭代中，步驟(二)所求的電壓值V_n達(dá)到預(yù)設(shè)的收斂誤差，此時(shí)計(jì)算得到的各節(jié)點(diǎn)的電壓值V_n即為所求電壓值，

S15：根據(jù)每一節(jié)點(diǎn)的電壓值繪制2D軸對(duì)稱非線性靜磁場中的磁勢(shì)云圖，如圖10所示為靜磁場中的磁勢(shì)云圖。

在本發(fā)明中，步驟S14中(二)對(duì)等效電路的求解可以使用節(jié)點(diǎn)電壓法進(jìn)行求解，其步驟為：

(一)計(jì)算得到矩陣YV＝I，其中Y為電路導(dǎo)納矩陣，由于迭代過程中，導(dǎo)納矩陣Y保持不變，僅需計(jì)算一次，V為待求節(jié)點(diǎn)電壓，I為節(jié)點(diǎn)電流；

(二)迭代第一步對(duì)矩陣Y進(jìn)行LU分解，即Y＝LU，其中L為單位下三角矩陣，U為上三角矩陣，之后的每一次迭代，無需計(jì)算此步，直接計(jì)算步驟(三)；

(三)使用公式V＝U^-1(L^-1I)求解節(jié)點(diǎn)電壓V。

在本發(fā)明中，步驟S14中(五)，入射過程中，每一非線性元件兩端電壓的求解是獨(dú)立的，在此使用多核并行技術(shù)同時(shí)對(duì)多個(gè)非線性元件兩端電壓進(jìn)行求解。

下面闡述本發(fā)明的有益技術(shù)效果：

相比傳統(tǒng)牛頓迭代法求解非線性有限元問題，使用本發(fā)明提供的傳輸線迭代法，能夠非常顯著的減少計(jì)算所用的時(shí)間。圖11～圖13對(duì)比了傳統(tǒng)牛頓迭代法與傳輸線迭代法的計(jì)算效率。圖11中，在單核計(jì)算下，傳統(tǒng)的牛頓迭代法計(jì)算時(shí)間幾乎是傳輸線迭代法的6倍，使用本發(fā)明提供的方法，求解速度得到了很大的提升；圖12中，隨著求解模型的分網(wǎng)單元增加，牛頓迭代法與傳輸線迭代法的求解時(shí)間比值不斷增大，說明出本發(fā)明能夠有效的處理有限元模型變復(fù)雜的情況；圖13顯示的是這兩種方法的單步求解時(shí)間對(duì)比，可見，本發(fā)明具有非常顯著的時(shí)間優(yōu)勢(shì)，能大幅度提高計(jì)算效率。

本發(fā)明提供的基于傳輸線迭代法的2D軸對(duì)稱非線性靜磁場模型的有限元求解方法解決了牛頓迭代法求解有限元非線性問題時(shí)帶來求解時(shí)間長，效率低的問題。本發(fā)明中的每一步驟均無需再次計(jì)算全局矩陣，只需要進(jìn)行一次全局矩陣的LU分解，然后即可重復(fù)使用，從而節(jié)省計(jì)算時(shí)間；同時(shí)，該方法非常適合使用并行算法進(jìn)行加速，能夠進(jìn)一步的加快有限元問題的求解。相對(duì)于傳統(tǒng)的牛頓迭代法，本發(fā)明在求解時(shí)間上有著很大的優(yōu)勢(shì)，有著廣闊的應(yīng)用前景。

本領(lǐng)域普通技術(shù)人員可以理解：附圖只是一個(gè)實(shí)施例的示意圖，附圖中的模塊或流程并不一定是實(shí)施本發(fā)明所必須的。

本領(lǐng)域普通技術(shù)人員可以理解：實(shí)施例中的裝置中的模塊可以按照實(shí)施例描述分布于實(shí)施例的裝置中，也可以進(jìn)行相應(yīng)變化位于不同于本實(shí)施例的一個(gè)或多個(gè)裝置中。上述實(shí)施例的模塊可以合并為一個(gè)模塊，也可以進(jìn)一步拆分成多個(gè)子模塊。

最后應(yīng)說明的是：以上實(shí)施例僅用以說明本發(fā)明的技術(shù)方案，而非對(duì)其限制；盡管參照前述實(shí)施例對(duì)本發(fā)明進(jìn)行了詳細(xì)的說明，本領(lǐng)域的普通技術(shù)人員應(yīng)當(dāng)理解：其依然可以對(duì)前述實(shí)施例所記載的技術(shù)方案進(jìn)行修改，或者對(duì)其中部分技術(shù)特征進(jìn)行等同替換；而這些修改或者替換，并不使相應(yīng)技術(shù)方案的本質(zhì)脫離本發(fā)明實(shí)施例技術(shù)方案的精神和范圍。