Skip to content

一体化构建方法

提示

构建方法基于张力结构找形找力一体化设计专利方法

以力密度为未知数的节点力平衡方程

建立结构模型,对结构节点、单元进行标号,构建张力结构关联矩阵Cs:

(1)Cs=[CCf]

其中:C(Rb×n)Cf(Rb×nf)分别为自由节点关联矩阵、固定节点关联矩阵,bnnf 分别为单元、自由节点、固定节点数量。假定单元 k 的两端点分别为节点 i、节点 j,则矩阵 Cskp 列为:

(2)cs(k,p)={1,p=i1,p=j0,pi,j

对每个自由节点,以单元力密度向量 q 为未知数,列力平衡方程:

(3){CTdiag(Csxs)q=pxCTdiag(Csys)q=pyCTdiag(Cszs)q=pz

其中:xs=[xxf]Tys=[yyf]Tzs=[zzf]Txyz(Rn) 分别为自由节点 xyz 坐标列向量:

(4){x=(x1,x2,,xn)Ty=(y1,y2,,yn)Tz=(z1,z2,,zn)T

xfyfzf(Rnf)分别为固定节点 xyz 坐标列向量:

(5){xf=(xn+1,xn+2,,xn+nf)Tyf=(yn+1,yn+2,,yn+nf)Tzf=(zn+1,zn+2,,zn+nf)T

pxpypz(Rn)分别为自由节点 xyz 方向外荷载:

(6){px=(px,1,xx,2,,xx,n)Tpy=(py,1,xy,2,,xy,n)Tpz=(pz,1,xz,2,,xz,n)T

q(Rb)=(q1,,qk,,qb)T,称为单元力密度列向量,元素 qk 为单元 k 的预应力 Nk 与长度 lk 的比值,diag() 表示以 () 为对角元素的方阵。

方程(3)进一步写为矩阵形式:

(7)Aq=p

其中:

(8)A(R3n×b)=[CTdiag(Csxs)CTdiag(Csys)CTdiag(Csys)]

称为平衡矩阵, p(R3n)=[pxpypz] ,称为节点外荷载列向量。

提示

现阶段sonew不考虑外荷载

单元预应力功能目标

单元预应力功能目标可表示为单元间力密度的相互关系:

(9)[α1,α2,,αb]q=0

如单元 uv 的力密度比值等于 α 可表示为:qu/qv=α,写为矩阵形式为:

(10)[0,,1,,α,,0]q=0

其预应力比值等于 α 可表示为:qu/qv=αlv/lu,写为矩阵形式为:

(11)[0,,1,,αlv/lu,,0]q=0

预应力功能目标写成矩阵形式:

(12)Ωq=ω

其中:Ω(Re×b) 为单元预应力功能目标系数矩阵,e 为单元预应力功能目标个数。

结构整体预应力水平

实际工程应用中,需要控制结构整体预应力水平,这里取单元力密度平方和为一常数,即:

(13)k=1bqk2=qc2

写为矩阵形式为:

(14)qTq=qc2

其中:qc2 为预应力总体水平系数。

提示

现阶段sonew只引入了此项预应力水平,且qc=1。这将在提交计算时,由程序自动施加。

以节点坐标为未知数的节点力平衡方程

对每个自由节点,以节点坐标 xyz 为未知数,列力平衡方程:

(15){CTdiag(q)Cx+CTdiag(q)Cfxf=pxCTdiag(q)Cy+CTdiag(q)Cfyf=pyCTdiag(q)Cz+CTdiag(q)Cfzf=pz

写成矩阵形式:

(16)[D000D000D][xyz]=[pxDfxfpyDfyfpzDfzf]

其中:D(Rn×n)=CTdiag(q)CDf(Rn×nf)=CTdiag(q)Cf,式(15)可进一步简化为:

(17)Dg[xyz]=[pxfpyfpxf]

其中:

(18)Dg(R3n×3n)=[D000D000D]

称为力密度矩阵,PxfPyfPzf 为考虑支座反力的节点外荷载向量。

节点坐标功能目标

节点坐标功能目标可表示为节点坐标的相互关系,其第 k 个约束关系可表示为:

(19)[ψk,1,ψk,2,,ψk,3n][xyz]=ψk

i 节点 x 方向坐标固定于 xi 可表示为:

(20)[0,,1,,0]x=xi

ij 节点 x 方向坐标相同可表示为:

(21)[0,,1,,1,,0]x=0

节点坐标功能目标写成矩阵形式:

(22)Ψ[xyz]=ψ

其中:Ψ(Rm×3n) 为节点坐标功能目标系数矩阵,m 为节点坐标功能目标个数,ω 为节点坐标功能目标相关系数。

基本方程组集

将式(7)(12)(14)(16)(22)组集,并注意到:diag(q)Cfxf=diag(Cfxf)q,可得下式:

(23)[CTdiag(q)C00CTdiag(Cfxf)0CTdiag(q)C0CTdiag(Cfyf)00CTdiag(q)CCTdiag(Cfzf)000CTdiag(Csxs)000CTdiag(Csys)000CTdiag(Cszs)Ψ0000Ω000qT][xyzq]=[pxpypzpxpypzψωqc2]

或:

(24)[D00Dfx0D0Dfy00DDfz000Ax000Ay000AzΨ0000Ω000qT][xyzq]=[pxpypzpxpypzψωqc2]

其中:

(25)D(Rn×n)=CTdiag(q)C(26){Dfx(Rn×b)=CTdiag(Cfxf)Dfy(Rn×b)=CTdiag(Cfyf)Dfz(Rn×b)=CTdiag(Cfzf)

ΨRm×3nΩRe×b

进一步简化得到找形找力一体化方程:

(27)Λζ=γ

其中:

(28)Λ(R(6n+m+e+1)×(3n+b))=[D00Dfx0D0Dfy00DDfz000Ax000Ay000AzΨ0000Ω000qT](29)ζ(R3n+b)=[xyzq](30)γ(R6n+m+e+1)=[pxpypzpxpypzψωqc2]

Λζγ 分别为找形找力一体化矩阵、坐标、系数

方程(27)将节点坐标、单元力密度两类完全不同的物理量统一到单一方程组,不再区分“形”与“力”,并综合节点坐标功能目标、单元预应力功能目标及结构整体预应力水平,形成形态构建一体化方程

方程的直接求解

方程(27)系数矩阵相对结构有限元刚度矩阵有以下特点:

  1. 行数大于列数。
  2. 为非对称矩阵。
  3. 一般情况下非正定。
  4. 含有未知数,方程(27)高度非线性。

由于初始结构并非在每个自由节点处均满足力平衡方程,方程(27)增广矩阵的秩大于系数矩阵的秩,即 rank([Λγ])>rank(Λ),方程(27)无解,只能求其近似解。常规做法采用最小二乘法进行迭代求解,其近似解为:

(31)ζ(k+1)=Λ(k)+γ

其中:Λ+(R(3n+b)×(6n+m+e+1))为矩阵 Λ 的广义逆,即

(32)Λ+=(ΛTΛ)1ΛT

直接求解存在两个问题:

  1. 方法对初始值非常敏感,迭代极易发散。
  2. 一些结构构建结果本身是多值的。例如平面索桁架,上下弦发生翻转后同样满足方程(27)。这使得迭代要么极易发散,要么收敛到错误结果。
平面索桁架上下弦翻转

平面索桁架上下弦翻转

方程的凸优化求解

为避免上述情况发生,必须改变求解策略。sonew将直接求解方程组问题转化为求解方程组误差最小的凸优化问题。同时附加不等式功能目标,进一步提高求解健壮性,并可避免多值问题。

提示

  1. 不等式约束将构建一道"铜墙铁壁",迭代值很难突破此界限。
  2. 通过不等式约束,将不需要的多值屏蔽,使可行集为凸集。可行集为凸集的凸函数,是一定有唯一极小值的。
  3. 简单模型、简单功能目标可不设置不等式约束。

以下简述凸优化方法。

对于非线性方程组:

(33)Ax=b

设:

(34)f(x)=AxbF=0.5fTf=0.5(Axb)T(Axb)

求解方程(33)可看作求函数 F 的极小值。当方程(33)有唯一解,F有唯一极小值,且极小值为0。

为实现不等式约束:

(35)Cx>d

提示

小于类型的不等式,可通过两边取负,转化为大于类型不等式。

采用内点罚函数法,引入罚函数:

(36)pi=1Cixdi

附带不等式约束的优化目标函数为:

(37)Q=F+αP

其中:

(38)P=0.5pTpp=[p1pipn]T

迭代产生最小值 Q(α) 及对应的 x(α),不断缩小 α 的值至零,即得到原问题的极小值 F 及对应的 x