test-posi.bdf

一个结构响应分析程序

$ NASTRAN input file created by the MSC MSC.Nastran input file
$ translator ( MSC.Patran 12.0.041 ) on December  23, 2005 at 13:51:33.
$ Direct Text Input for File Management Section
$ Normal Modes Analysis, Database
SOL 103
$ Direct Text Input for Executive Control
COMPILE SEMODES
ALTER 81           $'CALL.*PHASE1DR.*MAA'
$ 处理质量矩阵中的转动项
MATMOD MAA,,,,,/PARTNULL,MATRNULL/12/S,N,NONULL/3 $
add maa,matrnull/maa1//500.0
MATMOD MAA1,,,,,/AFILTER,/2///1/-1.0 $
add maa1,Afilter/negmaa/1.0/-1.0
add MAA1,negmaa/maa0/1.0/-1.0 $
$ 求出MAA矩阵的列宽NCOL
PARAML MAA0//'TRAILER'/2/S,N,NCOL $
TYPE PARM,,I,N,N2AA $
N2AA=NCOL*2 $
$ 声称单位阵I_I和零矩阵
MATGEN ,/I_I/1/ncol $
MATGEN ,/O_O/7/NCOL/NCOL $
$  形成装配矩阵AA
$  MK=-K/M
solve MAA0,KAA,,,/MK//-1/ $
$   生成组成AA的向量AACP
matgen ,/AACP/6/n2aa/ncol/ncol $
$   生成矩阵AA
merge O_O,MK,I_I,O_O,AACP,/AA $
$$  已经生成AA
$   生成矩阵BB  过程////在此设计载荷向量
solve MAA0,I_I,,,/MD  $使用全状态反馈
matgen ,/BBCP/6/n2aa/ncol/ncol $
matgen ,/BBCP1/6/ncol/ncol/ $
$  生成矩阵BB
merge O_O,I_I,,,BBCP1,BBCP/BB/ $
$   至此矩阵Ａ和Ｂ已经生成
$ 生成Q矩阵和R矩阵$$$$$$$$$$$$$$$
merge I_I,O_O,O_O,I_I,aacp,/Q_Q
add I_I,/R_R/1.00
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$ 求解riccati方程 A'K+KA+Q-KBR^B'K=0
$ 使用Newton迭代方法
$ 求解方法见叶庆凯 优化与最优控制中的计算方法 294页
$
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$ 1、计算belta的值，从而计算As
type parm,,rs,n,belta,rncol
type parm,,cs,n,beltac
norm MK/MKmax///s,n,maxmk
rncol=real(ncol)
belta=-0.5*rncol*maxmk
beltac=cmplx(belta)
MATGEN ,/I_Q/1/n2aa            $初始化为单位阵I_Q
add AA,I_Q/As/-1.0/beltac
$
$ 2、计算S=BR~Bt
trnsp BB/BBt
smpyad BB,R_R,BBt,,,/BBs/3 $
$
$ 3、求解方程的初值P0
type parm,,rs,n,errx=10.0
type parm,,i,n,ndo=0,mdo=0,kdo=0
add as,/as0
add as0,/as1
add bbs,/QQ/2.0
add qq,/qq0
add qq,/qq1
file as0=ovrwrt
file as1=ovrwrt
file qq0=ovrwrt
file qq1=ovrwrt

do while (ndo<1000 and errx>0.05e-35)
   SOLVE as0,,,,/INVas0/3
   add as0,invas0/as1/0.5/0.5
   trnsp invas0/invast
   smpyad invas0,qq0,invast,,,qq0/qqi/3
   add qqi,/qq1/0.5
   add As1,I_Q/as1Q
   norm as1q/maxas1///s,n,maxas1q
   errx=maxas1q
   copy as1/as0
   copy qq1/qq0
   ndo=ndo+1
   message //'maxas1 ='/errx
enddo
   add qq1,/pp0/0.5

trnsp as/ast
mpyad as,pp0,qq/aspp0
mpyad pp0,ast,aspp0/aspp0sub
norm aspp0sub/maxasp///s,n,maxaspp0
message //'the max aspp0 ='/maxaspp0
   MATMOD pp0,,,,,/P0NULL,p0RNULL/12/S,N,NONULL/3 $
   add pp0,p0rnull/p1//0.1e-5
   message //'ndo = '/ndo/ 'nonull='/nonull
   solve p1,I_Q,,,/kk0/$
   FILE KK0=OVRWRT
   FILE QQI=OVRWRT
   FILE ASPP0=OVRWRT
   FILE ASPP0SUB=OVRWRT
   FILE MAXASP=OVRWRT
   message //'ndo = '/ndo
$ 3、求解方程kk (8.50)
type parm,,rd,n,mrrorx=1000.0d0
type parm,,rd,n,kerrx0=10.0d100
type parm,,rd,n,kerrx1=10.0d0
type parm,,rd,n,krrx0=10.0d0
type parm,,rd,n,maxsub0=99.0d200
do while (mrrorx>0.5d-10 and mdo<1000)
   mpyad BBS,KK0,AA/aspt//-1
   message //'mdo= '/mdo
   trnsp aspt/asp
   smpyad kk0,bbs,kk0,,,Q_Q/qsp/3
   file QSP=OVRWRT
   FILE ASP=OVRWRT
   $求解asp*Pi+Pi*asp'=-qsp
   kdo=0
   kerrx0=10.0d1
   kerrx1=10.0d1
   krrx0=10.0d100
   do while (kerrx1>0.05d-30 and kdo<200) $16:59
      $SOLVE asp,,,,/INVasp/3
      DECOMP ASP/L,U,/ $
      FBS L,U,I_Q/invasp/ $
      add asp,invasp/asp1/0.5/0.5
      trnsp invasp/invaspt
      smpyad invasp,qsp,invaspt,,,qsp/qqi/3
      add qqi,/qq1/0.5
      add Asp1,I_Q/asp1Q
      FILE ASp1Q=OVRWRT
      norm asp1q/maxas1/////s,n,maxasp1q
      kerrx0=maxasp1q-kerrx1
      krrx0=kerrx0/kerrx1
      if (abs(krrx0)<0.1d-30) then $
         kdo=1000
      endif
      kerrx1=maxasp1q
      
      copy asp1/asp
      copy qq1/qsp
      kdo=kdo+1
      message //'maxasp1 ='/kerrx1
      file MAXAS1=OVRWRT
   enddo
      add qq1,/kk/0.5
      message //'kdo ='/kdo/'mdo = '/mdo
      add kk,kk0/subkk//-1.0
$验证
mpyad aa,kk,Q_Q/aspp0/1
mpyad kk,aa,aspp0/kaak
smpyad kk,BBS,kk,,,kaak/aspp0sub/3
norm aspp0sub/maxasp/////s,n,maxsub
$if (maxsub>maxsub0) then
$   mdo=10000
$endif
maxsub0=maxsub
file kaak=ovrwrt
      if (kdo>=900) then
         mrrorx=0.0d0
         message //'break out ! '
      else
         mrrorx=maxsub
         message //'maxsub = '/mrrorx
      endif
   copy kk/kk0
   mdo=mdo+1
enddo
message //'mdo = '/mdo
mpyad BB,kk0,/BK/1
matpch subkk,BK,kk0
$$$$$$$$验证riccati方程$$$$$$$
mpyad kk0,AA,Q_Q/kaq
mpyad AA,kk0,kaq/AKQ/1
smpyad kk0,bbs,kk0,,,/pbp/3
add akq,pbp/riccati/1.0/-1.0
norm riccati/ric///s,n,ricnum
message //'the riccati is = '/ricnum

CEND
SEALL = ALL
SUPER = ALL
TITLE = MSC.Nastran job created on 20-Dec-05 at 13:39:17
ECHO = NONE
$ Direct Text Input for Global Case Control Data
SUBCASE 1
$ Subcase name : Default
   SUBTITLE=Default
   METHOD = 1
   SPC = 2
   VECTOR(SORT1,REAL)=ALL
$ Direct Text Input for this Subcase
BEGIN BULK
PARAM    POST    0
PARAM    AUTOSPC YES
PARAM   PRTMAXIM YES
EIGRL    1       0.      200.    15      0                       MASS
$ Direct Text Input for Bulk Data
$ Elements and Element Properties for region : shell_qianlei
PSHELL   1       1       2.5     1               1
$ Pset: "shell_qianlei" will be imported as: "pshell.1"
CQUAD4   14602   1       4082    4220    4217    4118
CQUAD4   14605   1       4118    4217    4221    4116
CQUAD4   14606   1       4116    4221    4229    4097
CTRIA3   14609   1       4097    4229    4098
CQUAD4   14634   1       4013    4036    4208    4126
CQUAD4   14637   1       4126    4208    4212    4122
CQUAD4   14646   1       4014    4038    4219    4124
CQUAD4   14647   1       4124    4219    4223    4120
CQUAD4   14650   1       4120    4223    4227    4101
CQUAD4   14663   1       4018    4198    4252    4264
CQUAD4   14664   1       4264    4252    4251    4263
CQUAD4   14665   1       4263    4251    4203    4133
CQUAD4   14666   1       4023    4249    4201    4137
CTRIA3   14766   1       4101    4227    4102
CTRIA3   14770   1       4109    4207    4110
CTRIA3   14774   1       4129    4205    4130
CTRIA3   14775   1       4133    4203    4134
CTRIA3   14776   1       4137    4201    4138
CQUAD4   15476   1       4007    4188    4258    4268
CQUAD4   15477   1       4268    4258    4256    4267
CQUAD4   15478   1       4267    4256    4205    4129
CQUAD4   15479   1       4271    4210    4181    4073
CQUAD4   15480   1       4210    4271    4269    4214
CQUAD4   15481   1       4214    4269    4109    4207
CQUAD4   15494   1       4114    4112    4231    4233
CQUAD4   15495   1       4231    4112    4093    4234
CTRIA3   15497   1       4093    4234    4094
CQUAD4   15498   1       4244    4083    4085    4240
CQUAD4   15500   1       4114    4233    4237    4117
CQUAD4   15501   1       4117    4118    4217    4237