从ADINA中导出壳单元刚度矩阵,结果按照ADINA特定规范排列,通过MATLAB编写代码从ADINA导出的刚度矩阵文件读取并还原刚度矩阵
***GLOBAL STIFFNESS MATRIX, DIRECT MATRIX INPUT FILE*** | |||
---|---|---|---|
dmig-define stif-C:\Users\chenh\Desktop\ADINA\two_test_1 | |||
1 | X-T | 0.0000000000000E+00 | col |
1 | X-T | 0.1176271906564E+10 | |
1 | Y-T | 0.0000000000000E+00 | col |
1 | Y-T | 0.8191472696014E+09 | |
1 | X-T | -0.1785623184815E+09 | |
1 | Z-T | 0.0000000000000E+00 | col |
1 | Z-T | 0.8191472696014E+09 | |
1 | Y-T | -0.1190415456543E+09 | |
1 | X-T | 0.1785623184815E+09 | |
...... |
使用如下代码还原刚度矩阵
argv="plane.mtx"; % 文件名
infile=fopen(argv);
k=zeros(441*6); % 手动定义刚度矩阵大小
while ~feof(infile)
tline=fgetl(infile); % 按行读文件
a=textscan(tline,'%s');
if size(a{1},1)==4 % 确定当前刚度矩阵元素行数i
dof=a{1}(2);
dof=char(dof);
switch dof
case 'X-T'
dof=1;
case 'Y-T'
dof=2;
case 'Z-T'
dof=3;
case 'X-R'
dof=4;
case 'Y-R'
dof=5;
case 'Z-R'
dof=6;
end
indice=a{1}(1);indice=str2double(indice);
i=indice*6-6+dof;
elseif size(a{1},1)==3 % 确定当前刚度矩阵元素列数j
dof=a{1}(2);
dof=char(dof);
switch dof
case 'X-T'
dof=1;
case 'Y-T'
dof=2;
case 'Z-T'
dof=3;
case 'X-R'
dof=4;
case 'Y-R'
dof=5;
case 'Z-R'
dof=6;
end
indice=a{1}(1);indice=str2double(indice);
j=indice*6-6+dof;
k(i,j)=str2double(a{1}(3)); % 将当前刚度矩阵元素储存进刚度矩阵
end
end
k=tril(k,-1)+tril(k)'; % 下三角矩阵还原为整个刚度矩阵
运行程序,得到还原刚度矩阵如下:
1.176E9 | -1.179E8 | 1.179E8 | -3.670E-11 | -4.474E-10 | ... |
-1.179E8 | 8.192E8 | -1.190E8 | 5.834E7 | 8.251E7 | ... |
1.179E8 | -1.190E8 | 8.192E8 | 5.834E7 | 8.251E7 | ... |
-3.670E-11 | 5.834E7 | 5.834E7 | 2.918E7 | 2.104E3 | ... |
-4.474E-10 | 8.251E7 | 8.251E7 | 2.104E3 | 2.918E7 | ... |
... | ... | ... | ... | ... | ... |