OpenSees开发（二）源码分析——平面桁架静力有限元分析实例

本帖最后由 silyvin 于 2015-11-24 16:00 编辑

这是一个平面桁架静力分析算例，原始源代码位于 OpenSees2.3.0\EXAMPLES\Example1\main.cpp
因为字数限制，这里给一个下载，就不贴出来了，pan.baidu.com/s/1qWrdhyG

接下去一步一步解释代码：[code] // 创建一个有限元模型
Domain *theDomain = new Domain();[/code][code] // 创建4个节点，详细见说明1
Node *node1 = new Node(1, 2, 0.0, 0.0);
Node *node2 = new Node(2, 2, 144.0, 0.0);
Node *node3 = new Node(3, 2, 168.0, 0.0);
Node *node4 = new Node(4, 2, 72.0, 96.0); [/code]说明1：Node构造函数
位于OpenSees2.3.0\SRC\domain\node\Node.cpp
源码定义如下：
*****************************************************
(*Crd)(1) = Crd2;
。。。。。。
*****************************************************
参数tag为该节点的标签，指定给基类
: NomainComponent(tag,NOD_TAG_Node), NOD_TAG_Node是一个枚举值，表明该DomainComponent对象是一个节点类型；

ndof该节点的自由度，本例中，节点都为两个自由度；
Crd1, Crd2为该2维节点的坐标，赋于成员变量Crd，这是一个类数组的数据类型，创建了一个含该点坐标信息的数组。[code] // 将4个节点对象加入有限元模型中
// 如果两个node对象tag相同，则会返回失败
theDomain->addNode(node1);
theDomain->addNode(node2);
theDomain->addNode(node3);
theDomain->addNode(node4);[/code][code] // 创建一个弹性材料，见说明2
UniaxialMaterial *theMaterial = new ElasticMaterial(1, 3000);[/code]说明2：创建材料对象
*****************************************************
UniaxialMaterial *theMaterial = new ElasticMaterial(1,3000);
*****************************************************

UniaxialMaterial类官方说明：
opensees.berkeley.edu/OpenSees/api/doxygen2/html/classElasticMaterial.html
其中， ElasticMaterial 为 UniaxialMaterial 派生类

意为理想弹性材料
opensees.berkeley.edu/wiki/index.php/Elastic_Uniaxial_Material

构造函数
申明：
*****************************************************
ElasticMaterial(int tag, double E, double eta =0.0);
*****************************************************
实现：
*****************************************************
ElasticMaterial::ElasticMaterial(int tag, double e, doubleet)
:UniaxialMaterial(tag,MAT_TAG_ElasticMaterial),
trialStrain(0.0), trialStrainRate(0.0),
E(e), eta(et), parameterID(0)
{

}
*****************************************************
其中，第一个参数tag为标签，传递给基类构造函数，e为弹性模型，et为材料阻尼比，默认为0.[code] // 创建一个工况，编号为1，暂时未知
LoadPattern *theLoadPattern = new LoadPattern(1);
theDomain->addLoadPattern(theLoadPattern);

// 暂时未知这句话什么意思
theLoadPattern->setTimeSeries(new LinearSeries());

// 创建一个节点荷载向量
Vector theLoadValues(2);
theLoadValues(0) = 100.0;
theLoadValues(1) = -50.0;

// 第一个参数tag标签，第二个参数表明施加点荷载的节点tag，第三个参数是一个向量，表明在第一维度施加100个单位力，第二维度施加反方向50单位力
NodalLoad *theLoad = new NodalLoad(1, 4, theLoadValues);

// 将NodalLoad对象加入模型，1表示加入的工况编号
theDomain->addNodalLoad(theLoad, 1);

// 如果new NodalLoad后的节点编号未在模型中找到，返回失败
// 如果addNodalLoad第2个参数所表示的工况编号不存在，返回失败[/code]这里为了避免内存泄漏，也为了使代码的封装性更强，我更改了一部分代码：[code] AnalysisModel *theModel = new AnalysisModel();
EquiSolnAlgo *theSolnAlgo = new Linear();
StaticIntegrator *theIntegrator = new LoadControl(1.0, 1, 1.0, 1.0);
ConstraintHandler *theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8);
RCM *theRCM = new RCM();
DOF_Numberer *theNumberer = new DOF_Numberer(*theRCM);
BandSPDLinSolver *theSolver = new BandSPDLinLapackSolver();
LinearSOE *theSOE = new BandSPDLinSOE(*theSolver);

StaticAnalysis theAnalysis(*theDomain,
*theHandler,
*theNumberer,
*theModel,
*theSolnAlgo,
*theSOE,
*theIntegrator);[/code]改为：[code] // 分析对象封装
struct MyStaticAnalysis : public StaticAnalysis
{
ConstraintHandler *pConstraintHandler;
DOF_Numberer *pDOF_Numberer;
AnalysisModel *pAnalysisModel;
EquiSolnAlgo *pEquiSolnAlgo;
LinearSOE *pLinearSOE;
StaticIntegrator *pStaticIntegrator;

MyStaticAnalysis(Domain *theDomain) : StaticAnalysis(*theDomain,
*(pConstraintHandler = new PenaltyConstraintHandler(1.0e8,1.0e8)),
*(pDOF_Numberer = new DOF_Numberer(*(new RCM()))),
*(pAnalysisModel = new AnalysisModel()),
*(pEquiSolnAlgo = new Linear()),
*(pLinearSOE = new BandSPDLinSOE(*(new BandSPDLinLapackSolver()))),
*(pStaticIntegrator = new LoadControl(1.0, 1, 1.0, 1.0))) {}

~MyStaticAnalysis(){
delete pConstraintHandler;
delete pDOF_Numberer;
delete pAnalysisModel;
delete pEquiSolnAlgo;
delete pLinearSOE;
delete pStaticIntegrator;
}
};[/code][code] // 实例化分析模型对象
MyStaticAnalysis &theAnalysis = *(new MyStaticAnalysis(theDomain));

// perform the analysis & print out the results for the domain
int numSteps = 1;
theAnalysis.analyze(numSteps);

// 释放分析对象
delete &theAnalysis;

// 模型信息打印
opserr << *theDomain;

[/code][code] // 查看4节点两个自由度上的位移
Vector const &disp4node = theDomain->getNode(4)->getDisp();
printf("x4: %lf, z4: %lf\n", disp4node[0], disp4node[1]);

// 查看3个桁架单元的轴力
Information trussInfo;
for(int i=0; i<3; ++i)
{
Truss *pTruss = (Truss *)theDomain->getElement(i+1);
pTruss->getResponse(2, trussInfo);
printf("N%d: %lf\n", i+1, trussInfo.theDouble);
}

// Domain类的析构会释放加入domain的所有元素，所以node之类的对象不用单独析构
delete theDomain;[/code]编译——运行——屏幕输出：


第一自由度位移 0.530094，第二自由度位移-0.177894
杆件1轴力：43.94
杆件2轴力：-57.55
杆件3轴力：-55.31

与sap2000计算结果一致：





sap2000模型文件*.SDB(v14)和*.s2k文件，及修改后的源文件 first.cpp下载：
pan.baidu.com/s/1dDDKnb7

blog.csdn.net/silyvin/article/details/50008955
版权申明：本文为原创，仅发csdn和土木在线，禁止未经允许的转载。