彈性矩陣的案例
各向同性，橫觀各向同性，正交各向異性三種線彈性umat程序 ￥25
各向同性，橫觀各向同性，正交各向異性三種線彈性umat程序
1 各向同性
各向同性線彈性材料的彈性矩陣為：
式中拉梅常數的表達式為：
因此在編寫各向同性材料的umat時，需要兩個材料參數，在這里我們使用楊氏模量E和泊松比v。
2 橫觀各向同性
橫觀各向同性線彈性材料的彈性矩陣為：
并有關系式：
可見其彈性矩陣需要5個獨立的參數，為下列5個工程常數：
下標a代表軸向，下標t代表橫向。
3 正交各向異性
正交各向異性線彈性材料的彈性矩陣為：
并有關系式：
因此對于正交各向異性材料，其彈性矩陣需要9個工程常數來確定：
4 程序
使用Fortran90編寫umat程序。由于Abaqus默認的umat子程序為Fortran77，因此為了使用f90程序，使用命令：
abaqus make library=xxx.f90
該命令可以生成相應的后綴為obj的文件，之后使用該文件即可。使用上述方法可以避免使用Fortran77進行umat的編寫。
展開 ABAQUS用戶定義單元UEL與VUEL從入門到放棄系列2 ￥50
單元壓電耦合場的廣義本構方程表示如下：
即為：廣義應力=廣義彈性矩陣·廣義應變。其中廣義力中D為電位移，廣義應變中E為電場強度。
在廣義彈性矩陣中，C矩陣為力學場的彈性矩陣，e矩陣為壓電常數矩陣，右下角為介電常數矩陣。亦可展開如下表示（某種材料的參數，如果是特殊材料e矩陣中非0常數會更多或者更少，由材料本身決定）。
本構關系的張量表達式為：
其中，廣義應變的有限元格式可表示為：
其中，電場強度E為負的電勢梯度：
則廣義應變列陣記為：
那么，單元的剛度矩陣可以表示為：
其中Kuu為C3D8原本的剛度矩陣，Kuv與Kvu為壓電耦合剛度矩陣，Kvv為電場的廣義剛度矩陣。
ABAQUS中的C3D8和C3D8E都是做了一些剛度修正的，比如C3D8為了防止單元自鎖，采用了B-Bar方法，得出的剛度矩陣是介于C3D8和C3D8R之間的值，同樣的C3D8E也有一些類似的修正，以下我將提供一個不包含修正的版本，對ABAQUS剛度修正方法感興趣的朋友可以去拿去跟ABAQUS CAE對比。
展開 為什么復合材料層壓板設計中經常要求均衡性？
對于單向層合板，在材料坐標系下的應力應變關系如下，其中Q為彈性矩陣：
或
注意上述應力應變關系中Q16和Q26兩項均為0。
對于一個包含多個鋪層角度的 層壓板，其第k個鋪層，在層合板坐標系下單層的應力應變關系為：
其中， Q'為偏軸彈性矩陣，由 Q通過矩陣變換得到。
其展開形式可以表示為：
其中，
θ為鋪層角度。這里需要注意的是，偏軸彈性矩陣系數Q 11'、Q 22'、Q 66'、Q 12'四項均為 偶函數，+θ和-θ兩個鋪層對應的偏軸彈性系數相同，比如Q 11'(+θ)=Q 11'(-θ)。而Q 16'和Q 26'兩項是 奇函數，+θ和-θ兩個鋪層對應的這兩個偏軸彈性系數剛好 大小相等，符號相反，例如Q 16'(+θ)=-Q 16'(-θ)。
根據經典層合板理論，
先拋開B矩陣不管（以后再單獨說B矩陣的事）。
當一個層壓板中+θ和-θ鋪層數量相等時，即均衡鋪層時，A 16和A 26正負值成對出現，中性面上的正應變ε x0或ε y0在+θ和-θ層引起的剪力N xy大小相等、符號相反，互相抵消。也就是面內的正應變不會引起附加的剪切變形， 不存在拉剪耦合效應。
作用在層壓板上的力
反之，如果一個層壓板中+θ和-θ鋪層數量不相等時，面內正應變引起的剪力無從抵消，層板在發生拉伸/壓縮變形的同時，必然還有附加的剪切變形，即存在 拉剪耦合效應，如下圖所示。同樣滴，當層壓板中面有剪切變形γ xy時，除了引起剪力N xy之外，還會引起軸向力N x、N y，即存在 剪拉耦合效應。
展開 彈性力學基本方程的矩陣形式
平衡方程
彈性體V域內任一點沿坐標軸x,y,z方向平衡方程
其中 ， ， 為單位體積的體積力在x,y,z方向的分量
平衡方程矩陣形式
其中
是體積向量，
【姊妹篇】為什么復合材料層壓板設計中經常要求對稱性？
上一節曾經描述過，對于單向層合板，在材料坐標系下的應力應變關系如下，其中Q為彈性矩陣：
或
對于一個包含多個鋪層角度的
層壓板
，其第k個鋪層，在層合板坐標系下單層的應力應變關系為：
其中，
Q'為偏軸彈性矩陣，由
Q通過矩陣變換得到。
其展開形式可以表示為：
其中，
θ為鋪層角度。這里需要注意的是，偏軸彈性矩陣系數Q
11'、Q
22'、Q
66'、Q
12'四項均為
偶函數，+θ和-θ兩個鋪層對應的偏軸彈性系數相同，比如Q
11'(+θ)=Q
11'(-θ)。而Q
16'和Q
26'兩項是
奇函數，+θ和-θ兩個鋪層對應的這兩個偏軸彈性系數剛好
大小相等，符號相反，例如Q
16'(+θ)=-Q
16'(-θ)。
根據經典層合板理論，
在講均衡性的時候重點看的是A16、A26、D16、D26。本期重點看一下耦合剛度矩陣B.
根據上式中B矩陣的表達式，當相同材料、相同鋪層角度的兩個鋪層對稱分布在層壓板中性面兩側時，相同材料保證了Q陣的一致，相同鋪層角度保證了Q'矩陣的一致。假定層壓板中一共有N個鋪層，第k鋪層與第N-k層對稱，如下圖所示。兩個對稱鋪層，無論角度是多少，其累加之后的B矩陣都是0。
圖2 層板厚度方向幾何分層示意圖
圖3 B矩陣為0
B矩陣為0，表示面內的拉伸、壓縮、剪切變形不會誘導出面外的彎曲、扭轉變形，即不存在面內-面外、面外-面內耦合。
展開 彈性力學基本方程之應力應變的矩陣形式
應力列陣（應力向量）
應力分量正負號規定
如果某一個面的外法線方向與坐標軸的正方向一致，這個面上的應力分量就以沿坐標軸正方向為正，與坐標軸反方向為負；相反如果某一個面的外法向方向與坐標軸的負方向一致，這個面上的應力分量就以沿坐標軸負方向為正，與坐標軸同向為負。應力分量及正方向如圖1
Abaqus 中常用的復合材料本構介紹【圖文+視頻】
正交各向異性Orthotropic
這類材料本構在CAE中需要輸入的是彈性矩陣系數，計算公式如下：
完全各向異性材料Anisotropic
三維完全各向異性材料：直接指定彈性矩陣，無對稱面，共21個獨立彈性系數。
點擊播放視頻
該視頻主要講解了Abaqus 中常用的復合材料本構。
回顧經典層合板理論
另外，基于Kirchhoff假設，變形后的法線仍然垂直于中面，變形后法線仍為直線且長度不變，不考慮沿厚度方向的剪切變形，則有：
3
應力應變關系
對于單向層合板，在材料坐標系下的應力應變關系為，
寫成柔度矩陣形式為，
對柔度矩陣求逆可以得到彈性矩陣Q，彈性矩陣及柔度矩陣中的系數表達式可查閱任意一本復合材料力學書。
Ansys材料參數的定義問題
相關命令,例如:
mp,ex,1,300e9
mp,ey,1,200e9
mp,nuxy,1,0.25
mp,gxy,1,170e9
…
3.anisotropic:各向異性材料:
各向異性材料定義起來較為復雜,這里我只作些簡單的說明,更詳細的資料,大家可以去看幫助.對于各向異性彈性材料的定義,需要定義彈性系數矩陣,這個矩陣是一個對稱正定陣,因而輸入的值一定要為正值.
彈性常數矩陣如下圖所示,各向異性體只有21個獨立的彈性常數,因而我們也就只需輸入21個參數即可,而且對于二維問題,彈性常數縮減為10個.彈性系數矩陣可以用剛度或柔度兩種形式來定義,自己根據情況選用,輸入的時候,可以通過菜單或者TB命令的TBOPT選項來控制.
相關的命令流,例如:
tb,anel,1
tbdata,1, 110e6, 120e6, 130e6, 140e6, 150e6, 160e6
tbdata,7, 220e6, 230e6, 240e6, 250e6, 260e6
tbdata,12, 330e6, 340e6, 350e6, 360e6
tbdata,16, 440e6, 450e6, 460e6
tbdata,19, 550e6, 560e6
tbdata,21, 660e6
另:需注意一下各個參數的編號順序和起始位置,不要搞錯了,輸入的時候,是按照上三角陣來錄入的,即:D11,D12,D13,D14,D15,D16,D22,D23…
展開 有限元單元矩陣計算流程
集合單元矩陣流程 (Assem.for)
◎ 確定各計算參數，如質量密度為sp, 調用子程序element_vD_plane, 根據問題類型mprob, 確定單元彈性矩陣vD
◎ 單元質量矩陣和單元剛度矩陣清零。
◎ 在單元內部，按積分點循環，計算各積分點上的矩陣。
◎ 將各積分點上的矩陣，按積分點的權重進行加權平均，即可得到各單元矩陣。
ANSYS APDL中的壓電分析
設應力為T，加于截面積A的壓電陶瓷片上，其所產生的應變為S，則根據胡克定律，應力T與應變S之間有如下關系
S=sT ，
T=cS
式中，S為彈性順度常數，單位為m2/N；C為彈性勁度常數，單位為N/m2
但是，任何材料都是三維的，即當施加應力于長度方向時，不僅在長度方向產生應變，寬度與厚度方向上也產生應變。設有如圖1-2所示的薄長片，其長度沿1方向，寬度沿2方向。沿1方向施加應力T1，使薄片在1方向產生應變S1，而在方向2上產生應變S2，由（1-5）式不難得出
S1=S11T1 ，
S2=S12T1
上面兩式彈性順度常數S11和S12之比，稱為泊松比，它表示橫向相對收縮與縱向相對伸長之比。
由于壓電陶瓷存在壓電效應，因此，壓電陶瓷樣品在不同的電學條件下具有不同的彈性順度常數。在外電路的電阻很小相當于短路，或電場強度E=0的條件下測得的稱為短路彈性順度常數，記作SE。在外電路的電阻很大相當于開路，或電位移D=0的條件下測得的稱為開路彈性順度常數，記作SD。由于壓電陶瓷為各向異相性體，因此共有下列10個彈性順度常數：
SE11，SE12，SE13，SE33，SE44，SD11，SD12，SD13，SD33，SD44。
同理，彈性勁度常數也有10個：
CE11，CE12，CE13，CE33，CE44，CD11，CD12，CD13，CD33，CD44
ANSYS中，彈性系數矩陣為6×6矩陣（對2-D模型是4×4矩陣），它說明剛度系數（[c]矩陣）或柔度系數([s]矩陣)。
展開 
各向同性彈塑性本構的vumat源代碼：通過修改umat ￥99
彈性矩陣
real(8) :: E,nu
real(8) :: De(6,6)
! 計算全應變
strain_n1 = strain_n0 + strain_inc
! 計算彈性矩陣
E = 9.0 * kappa * mu / (3.0 * kappa + mu)
nu = ( 3.0 * kappa - 2.0 * mu ) / (3.0 * kappa + mu) / 2.0
De = plastic_iso_compute_elasD(E,nu)
! 讀取n時刻的狀態變量
alpha_n0 = statev_n0(1)
strain_pl_n0 = statev_n0(2:7)
Dp_0 = statev_n0(8)
!
展開 單元剛度矩陣
由（27a）式定義的單元剛度矩陣，由于應變B對于3結點三角形單元是常量陣，因此有
代入彈性矩陣D和應變矩陣B后，它的任一分塊矩陣可表示成
由(34)式立即可以得到
由此可見單元剛度矩陣是對稱矩陣。
為了進一步理解單元剛度矩陣的物理意義，我們同樣可以利用最小位能原理建立一
個單元的平衡，這就得到
Pe是單元結點載荷，當然應當包括其它相鄰單元對該單元的作用力?，F把ae、Pe順序表示為
這是單元結點平衡方程，每個結點在x和y方向上各有一個平衡方程，3個結點共有六個平衡方程。方程左端是通過單元結點位移表示的單元結點內力，方程右端是單元結點外載。
令a1=1(ui=1),a2=a3=…＝a6=0
由（38）式可以得到
（39）式表明，單元剛度矩陣第一列元素的物理意義是：a1＝1。其他結點位移都為零時，需要在單元各結點位移方向上施加結點力的大小。當然，單元在這些結點力作用下處于平衡，因此在x和y 方向上結點力之和為零：
對于單元剛度矩陣中其他列的元素也可用同樣的方法得到它們的物理解釋。因此單元剛度矩陣中任一元素Kij物理意義為：當單元的第j個結點位移為單位位移而其他結點位移為零時，需在單元第j個結點位移方向上施加的結點力的大小。單元剛性大，則使結點產生單位位移所需施加的結點力就大。因此單元剛度矩陣中的每個元素反映了單元剛性的大小，稱為剛度系數。
展開 基于Runge-Kutta算法的硬化土模型二次開發
2 應力更新計算方法優化
ABAQUS中編寫自定義本構模型需要推導出在笛卡爾三維（或二維）坐標系下單元積分點的雅可比矩陣DDSDDE，即求解出該積分點各應力增量對應變增量的變化率?σ/?ε，并據此對該增量步的總應力和狀態變量進行更新、輸出，通過接口返回ABAQUS主程序[6]。
2.1 變形階段劃分與試探應力
對任意傳入的應變增量，先假定其全部為彈性應變，由廣義胡克定律可得材料的彈性剛度矩陣：
式中：[De]為彈性剛度矩陣；λ為拉梅常數；G為剪切模量。
式中：E為彈性模量；μ為泊松比。
式中：{?σn+1}為第n+1步的應力增量張量；{?σen+1}為第n+1步的彈性應力增量張量。
式中：為試探應力；為第n步的應力張量；{?σn+1}為第n+1步的應力張量。
將代入式（1）和式（9），若，則增量步全程處于彈性階段；若，則增量步從彈性階段過渡至彈塑性階段；若，則增量步全程處于彈塑性階段。
2.2 應力更新計算
在彈塑性階段，材料在從到的過程中必然會經過屈服應力狀態，在從到此段內只發生彈性變形，僅須用彈性剛度矩陣[De]更新應力即可。在到這段則發生彈塑性變形，須求得彈塑性剛度矩陣[Dep]進行求解。因此，為求解整個增量步的應力變化，須確定在屈服面上的位置，可設：
式中：r為比例因子；為第n+1步的屈服應力張量。
則：
由在屈服面上可得：
式中：為關于判斷應力的函數；F (σn+r·?σn+1)為關于應力及參數r的函數；F(r)為關于r的函數。
展開 關于HuangUMAT代碼中變量的中文翻譯
需要更多參數來表征單晶本構規律的用戶，例如Zarka提出的框架應該使NSTATV大于（或等于）這些參數的數量NPARMT加上滑動系統總數的九倍，NSLPTL，再加上五
NSTATV >= NPARMT+9*NSLPTL+5
(2) 如果考慮潛在硬化，切線剛度矩陣通常不是對稱的。用戶必須在輸入文件中的*USER MATERIAL卡上聲明“UNSYMM”。
參數ND確定此子程序中數組的維數。當前選擇150是立方晶體的上限，最多可激活三組滑移系統。用戶可以將參數ND減少到任意數量，只要大于或等于所有組中的滑移系統總數。 例如，如果{110}<111>是唯一可能激活的滑移系集合，則ND可被視為十二(12)。
NSLIP——每套內滑移系數量。
SLPDIR——滑移方向（初始狀態的單位向量）。
SLPNOR——滑移面的法線（初始狀態下的單位法線）。
SLPDEF——滑移變形張量（施密特因子）。
SLPSPN——滑移自旋張量（僅在有限旋轉時需要）。
DSPDIR——滑移方向的增量。
DSPNOR——滑移平面法線的增量。
DLOCAL——局部立方晶系中的彈性矩陣。
D——全局系統中的彈性矩陣。
ROTD——旋轉矩陣將DLOCAL轉換為D。
ROTATE——旋轉矩陣，全局系統中立方晶體的[100]、[010]和[001]的方向余弦。
FSLIP——滑移系的剪切應變率。
DFDXSP——FSLIP對x=TAUSLP/GSLIP的導數，其中TAUSLP是分剪切應力，GSLIP是當前強度。
DDEMSD——彈性模量張量與滑移變形張量的雙點積，僅對于有限旋轉，滑移自旋張量與應力的點積。
H——自硬化和潛在硬化矩陣。
DDGDDE——滑移系中剪切應變增量對應變增量的導數。
展開 