OptiStruct输出结构的刚度与质量矩阵

资料

1. https://help.altair.com/hwsolvers/os/topics/solvers/os/param_post_bulk_r.htm
2. https://help.altair.com/hwsolvers/os/topics/solvers/os/output_io_r.htm
3. https://help.altair.com/hwsolvers/os/topics/solvers/os/full_mat_file_r.htm

概述

已知，OptiStruct 有三种方法可以输出结构的刚度矩阵与质量矩阵。

方法 1：使用下面的参数可以生成稀疏格式的或者稠密格式的刚度矩阵与质量矩阵，结果存储在.full.mat文件，这是一个 ASCII 文件

$ 只能选一种
OUTPUT,MATRIX,FULL,SPARSE
OUTPUT,MATRIX,FULL,DENSE

方法 2：使用参数PARAM, POST, -5，这样会生成.k.op2与 .m.op2 文件，里面含有刚度矩阵与质量矩阵信息，但是目前暂无读取该 op2 文件的方法。

方法 3：使用超单元方法。

方法 1 OUTPUT

获得.full.mat 文件

1. 在 HyperMesh 中完成特征值分析设置，附件是一个特征值分析的 fem，设置方法可参考 https://19970622.xyz/posts/57316/
2. 在 fem 第 1 行添加OUTPUT,MATRIX,FULL,SPARSE
3. 使用 OptiStruct 求解，获得.full.mat文件
4. 需要指出 OptiStruct 不会对矩阵进行半带宽优化，文件也给出了刚度矩阵每个自由度含义；如果设置约束，导出的矩阵将是划行划列后的

读取.full.mat 文件

文件格式见：https://help.altair.com/hwsolvers/os/topics/solvers/os/full_mat_file_r.htm

读取函数：Read_KM_From_Optistruct.m

优化后的函数：Read_KM_From_Optistruct_Fast.m

优化后的函数：3993DOFs 需要 1.5 秒，336633DOFs 需要 143 秒

优化方向是使用sparse(i,j,v)

下面给出 Read_KM_From_Optistruct.m，注意waitbar会耗费一定时间

function Model = Read_KM_From_Optistruct(filName)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ref: https://www.yuque.com/xdd1997/ek3kug/xmmw08php2yi79yr                   %
% PWD: cawk                                                                    %
% update:2024-11-11                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(filName,'r');
TEXT = textscan(fid,'%s','Delimiter','\n','Whitespace','');
TEXT=char(TEXT{1});
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii = 1:size(TEXT,1)
    str00 = TEXT(ii,1:4);
    if     strcmpi(str00,'STIF')
        K_Ben = ii;
    elseif strcmpi(str00,'MASS')
        M_Ben = ii;
        K_End = ii-2;
    elseif strcmpi(str00,'DOFS')
        DOF_Ben = ii;
        if exist('M_Ben','var')
            M_End = ii-2;
        else
            K_End = ii-2;
        end
    end
    DOF_End = ii-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K_TEXT = TEXT(K_Ben:K_End,:);
DOF_TEXT = TEXT(DOF_Ben:DOF_End,:);
if exist('M_Ben','var')
    M_TEXT = TEXT(M_Ben:M_End,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% K
fprintf(string(datetime('now'))+":开始读刚度矩阵\n")
tmp = textscan(K_TEXT(1,:),'%s%d%d%d%d%d%s');
NCOL = tmp{5};
NROW = tmp{6};
if tmp{4}==1
    K = zeros(NROW,NCOL);
elseif tmp{4} ==3
    K = sparse(NROW,NCOL);
end
TEXT_2 = K_TEXT(2:end,:);
oneRow = reshape(TEXT_2',1,[]);
NCharOneRow = size(TEXT_2,2);
NOneRowText = length(oneRow);
POS = 0;

h = waitbar(0);
while true
    ICOL = str2double(oneRow(POS+[1:8]));
    rowRange = str2double(oneRow(POS+[9:16])):str2double(oneRow(POS+[17:24]));
    POS = POS + NCharOneRow;
    ValueText = oneRow([1:length(rowRange)*16]+POS);
    K(rowRange,ICOL) = double(string(reshape(ValueText,16,[])'));
    POS = POS+ceil(length(rowRange)/8)*NCharOneRow;

    if POS == NOneRowText
        delete(h);
        break
    else
        waitbar(ICOL/NCOL, h, sprintf('读刚度矩阵进度：%d/%d', ICOL,NCOL));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% K


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M
if exist('M_Ben','var')
    fprintf(char(string(datetime('now'))+":开始读质量矩阵\n"))
    tmp = textscan(M_TEXT(1,:),'%s%d%d%d%d%d%s');
    NCOL = tmp{5};
    NROW = tmp{6};
    if tmp{4}==1
        M = zeros(NROW,NCOL);
    elseif tmp{4} ==3
        M = sparse(NROW,NCOL);
    end
    TEXT_2 = M_TEXT(2:end,:);
    oneRow = reshape(TEXT_2',1,[]);
    NCharOneRow = size(TEXT_2,2);
    NOneRowText = length(oneRow);
    POS = 0;

    h = waitbar(0);
    while true
        ICOL = str2double(oneRow(POS+[1:8]));
        rowRange = str2double(oneRow(POS+[9:16])):str2double(oneRow(POS+[17:24]));
        POS = POS + NCharOneRow;
        ValueText = oneRow([1:length(rowRange)*16]+POS);
        M(rowRange,ICOL) = double(string(reshape(ValueText,16,[])'));
        POS = POS+ceil(length(rowRange)/8)*NCharOneRow;

        if POS == NOneRowText
            delete(h)
            break
        else
            waitbar(ICOL/NCOL, h, sprintf('读质量矩阵进度：%d/%d', ICOL,NCOL));
        end
    end
else
    M = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DOF
fprintf(char(string(datetime('now'))+":开始读自由度\n"))
tmp = textscan(DOF_TEXT(1,:),'%s%d%d%d%d%d%s');
NCOL = tmp{5};
NROW = tmp{6};
if tmp{4}==1
    DOF = zeros(NROW,NCOL);
elseif tmp{4} ==3
    DOF = sparse(NROW,NCOL);
end
TEXT_2 = DOF_TEXT(2:end,1:64);
oneRow = reshape(TEXT_2',1,[]);
NCharOneRow = size(TEXT_2,2);
NOneRowText = length(oneRow);
POS = 0;

while true
    ICOL = str2double(oneRow(POS+[1:8]));
    rowRange = str2double(oneRow(POS+[9:16])):str2double(oneRow(POS+[17:24]));
    POS = POS + NCharOneRow;
    ValueText = oneRow([1:length(rowRange)*8]+POS);
    DOF(rowRange,ICOL) = double(string(reshape(ValueText,8,[])'));
    POS = POS+ceil(length(rowRange)/8)*NCharOneRow;

    if POS == NOneRowText
        break
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DOF

Model.K = K;
Model.M = M;
Model.DOF = DOF;

方法 2 POST

获得 op2 文件

1. 在 HyperMesh 中完成特征值分析设置
2. 在 fem 文件的BEGIN BULK后面添加PARAM,POST,-5
3. 使用 OptiStruct 求解，获得.k.op2与 .m.op2 文件

读取 op2 文件

暂时没有办法。

pyNastran最新版本v1.4.1也不能读取这个 op2 文件，输出如下：

方法 3 超单元法

1. 划分网格—设置材料—设置单元属性
2. 右键创建 Load Collectors—Analysis—constraints—node—all—load types=ASET—creat
3. 右键创建 Load Step—Analysis type=Generic，注意不要配置 SPC
4. Analysis—OptiStruct—export options=all—run options=analysis—OptiStruct
5. 刚度矩阵在文件*_AX.pch中

附：OptiStruct 特征值分析 fem

OUTPUT,MATRIX,FULL,SPARSE
DISPLACEMENT(OUTPUT2,REAL) = YES

SUBCASE        1
ANALYSIS MODES
  METHOD(STRUCTURE) =        1

BEGIN BULK
PARAM,POST,-5

GRID           1        -1.0    -1.0    -1.0
GRID           2        1.0     -1.0    -1.0
GRID           3        1.0     1.0     -1.0
GRID           4        -1.0    1.0     -1.0
GRID           5        -1.0    -1.0    1.0
GRID           6        1.0     -1.0    1.0
GRID           7        1.0     1.0     1.0
GRID           8        -1.0    1.0     2.0

CHEXA          1       1       1       2       3       4       5       6
+              7       8

PSOLID         1       1
MAT1           1210000.0        0.3     7.85-9
EIGRL          1                      10
ENDDATA