PlaneTruss
。
m
function
PlaneTruss
%本程式用於求解平面桁架問題的節點力和節點位移
%使用前請將本檔案（‘PlaneTruss。m’）與‘myTruss。m’放置在同一資料夾中
%編寫時間：2020年5月25日-28日
setData
；
%初始化資料
SovelModel
；
%求解
getRes
；
%輸出結果
return
%Step1
function
setData
%讀取結構引數
global
Node
Element
Material
BC1
F
［
Node
，
Element
，
Material
，
BC1
，
F
］=
feval
（
‘myTruss’
）；
return
%Step2
function
SovelModel
global
Node
Element
nodeNum
eleNum
F
totalK
Delta
Reaction
axisF
；
［
nodeNum
，
~
］=
size
（
Node
）；
［
eleNum
，
~
］=
size
（
Element
）；
totalK
=
sparse
（
nodeNum
*
2
，
nodeNum
*
2
）；
for
i
=
1
：
1
：
eleNum
［
localKe
，
c
，
s
］=
gEStiffness
（
i
）；
Ke
=
gTransK
（
localKe
，
c
，
s
）；
gTStiffness
（
i
，
Ke
）；
end
modK
=
gBoundary
；
Delta
=
modK
\
F
；
Reaction
=
totalK
*
Delta
；
axisF
=
gElementF
；
return
%Step3
function
getRes
%輸出結果
global
nodeNum
eleNum
Delta
Reaction
axisF
fprintf
（
‘ 節點位移\n’
）；
fprintf
（
‘節點號 x方向位移 y方向位移\n’
）；
for
i
=
1
：
1
：
nodeNum
fprintf
（
‘%3d %10。6f %10。6f\n\n’
，
i
，
Delta
（（
i
-
1
）
*
2
+
1
，
1
）
+
0
，
Delta
（（
i
-
1
）
*
2
+
2
，
1
）
+
0
）；
end
fprintf
（
‘ 節點外力\n’
）；
fprintf
（
‘節點號 x方向外力 y方向外力\n’
）；
for
j
=
1
：
1
：
nodeNum
fprintf
（
‘%3d %10。6f %10。6f\n’
，
j
，
Reaction
（（
j
-
1
）
*
2
+
1
，
1
）
+
0
，
Reaction
（（
j
-
1
）
*
2
+
2
，
1
）
+
0
）；
end
fprintf
（
‘ 杆件軸力\n’
）；
fprintf
（
‘單元號 軸力\n’
）；
for
k
=
1
：
1
：
eleNum
fprintf
（
‘%3d %10。6f\n’
，
k
，
axisF
（
k
，
1
）
+
0
）；
end
return
function
［localKe，c，s］
=
gEStiffness
（
i
）
%計算區域性座標系下的單元剛度陣
global
Material
Node
Element
E
=
Material
（
Element
（
i
，
4
），
1
）；
% I=Material（Element（i，4），2）； %桁架暫時用不到I
A
=
Material
（
Element
（
i
，
4
），
3
）；
x1
=
Node
（
Element
（
i
，
2
），
2
）；
y1
=
Node
（
Element
（
i
，
2
），
3
）；
x2
=
Node
（
Element
（
i
，
3
），
2
）；
y2
=
Node
（
Element
（
i
，
3
），
3
）；
L
=
sqrt
（（
x2
-
x1
）
^
2
+
（
y2
-
y1
）
^
2
）；
c
=
（
x2
-
x1
）
/
L
；
s
=
（
y2
-
y1
）
/
L
；
localKe
=
（
E
*
A
/
L
）
。*
［
1
0
-
1
0
0
0
0
0
-
1
0
1
0
0
0
0
0
］
；
return
function
Ke
=
gTransK
（
localKe，c，s
）
%轉換為整體座標系下的單元矩陣
global
T
T
=［
c
s
0
0
-
s
c
0
0
0
0
c
s
0
0
-
s
c
］
；
Ke
=
T
。
‘*
localKe
*
T
；
return
function
gTStiffness
（
ie，Ke
）
%整合單元剛度陣到整體剛度陣
global
Element
totalK
for
i
=
1
：
1
：
2
for
j
=
1
：
1
：
2
for
m
=
1
：
1
：
2
for
n
=
1
：
1
：
2
x
=（
i
-
1
）
*
2
+
m
；
y
=（
j
-
1
）
*
2
+
n
；
X
=（
Element
（
ie
，
i
+
1
）
-
1
）
*
2
+
m
；
Y
=（
Element
（
ie
，
j
+
1
）
-
1
）
*
2
+
n
；
totalK
（
X
，
Y
）=
totalK
（
X
，
Y
）
+
Ke
（
x
，
y
）；
end
end
end
end
return
function
modK
=
gBoundary
（）
%處理邊界條件
%採用乘大數法
global
BC1
totalK
F
modK
=
totalK
；
［
i
，
~
］=
size
（
BC1
）；
for
x
=
1
：
1
：
i
m
=（
BC1
（
x
，
1
）
-
1
）
*
2
+
BC1
（
x
，
2
）；
F
（
m
，
1
）=
F
（
m
，
1
）
*
BC1
（
x
，
3
）
*
10e30
；
modK
（
m
，
m
）=
modK
（
m
，
m
）
*
1e30
；
end
return
function
axisF
=
gElementF
（）
%求解杆單元軸力
global
Element
Delta
T
［
eleNum
，
~
］=
size
（
Element
）；
axisF
=
sparse
（
eleNum
，
1
）；
delta
=
zeros
（
4
，
1
）；
for
i
=
1
：
1
：
eleNum
［
localKe
，
c
，
s
］=
gEStiffness
（
i
）；
T
=［
c
s
0
0
-
s
c
0
0
0
0
c
s
0
0
-
s
c
］
；
delta
（
1
，
1
）=
Delta
（
2
*
Element
（
i
，
2
）
-
1
，
1
）；
delta
（
2
，
1
）=
Delta
（
2
*
Element
（
i
，
2
），
1
）；
delta
（
3
，
1
）=
Delta
（
2
*
Element
（
i
，
3
）
-
1
，
1
）；
delta
（
4
，
1
）=
Delta
（
2
*
Element
（
i
，
3
），
1
）；
delta
=
T
*
delta
；
Force
=
localKe
*
delta
；
%讀取區域性座標系下的軸力
%結果是雙零一正一負的單列矩陣，注意力的方向
for
x
=
1
：
1
：
4
if
（
Force
（
x
）
~=
0
）
axisF
（
i
）=
-
Force
（
x
）；
break
；
end
end
end
return
myTruss
。
m
function
［Node，Element，Material，BC1，F］
=
myTruss
（）
%整體座標系取水平為x方向，豎直為y方向的笛卡爾直角座標系
%節點號/x座標/y座標
Node
=［
1
0
1
2
1
1
3
0
0
4
1
0
］；
%單元編號/節點號/節點號/材料種類
Element
=［
1
1
2
1
2
1
3
1
3
2
4
1
4
1
4
1
5
3
4
1
］；
%材料彈模/截面慣性矩/面積
Material
=［
1
1
1
］；
%第一類邊界條件，直接給出β值
%節點號/自由度種類（1水平2豎直）/β值（位移限定值）
BC1
=［
3
1
0
3
2
0
4
1
0
4
2
0
］；
%輸入已知的節點外力
［
nodeNum
，
~
］=
size
（
Node
）；
F
=
sparse
（
nodeNum
*
2
，
1
）；
F
（
3
，
1
）=
2
；
F
（
4
，
1
）=
-
1
；
end
