#2つの片持ち梁の弾塑性解析(STA10)
# どちらも応力-ひずみ曲線は、多直線近似しているものの、
# 一方は弾性体、もう片方は弾塑性体なので、除荷後の挙動に違いが見られる。

DEBUT();

#(1)応力-ひずみ曲線の多直線近似((ひずみ,応力)の並び)
elast_pl=DEFI_FONCTION(NOM_PARA='EPSI',NOM_RESU='SIGM',VALE=(0.001,1.0, #初期降伏点
                             0.1,1.99,
                             1,10.99,
                             2,20.99,
                             ),);

#(2)読み込みメッシュファイルの形式
MAIL=LIRE_MAILLAGE(FORMAT='MED',);

#(3)材料定数
E1=DEFI_MATERIAU(ELAS=_F(E=1000,                # MPa :ヤング率
                         NU=0.499,),            #ポアソン比
                 TRACTION=_F(SIGM=elast_pl,),); #(1)の名前

#(4)全てのものが3次元の問題ですよと
MODE=AFFE_MODELE(MAILLAGE=MAIL,
                 AFFE=_F(TOUT='OUI',
                         PHENOMENE='MECANIQUE',
                         MODELISATION='3D',),);

#(5)要素グループごとに与える材料定数
MATE=AFFE_MATERIAU(MAILLAGE=MAIL,
                   AFFE=_F(TOUT='OUI',   #すべての要素グループに適用
                           MATER=E1,),); #(3)の名前

#(6)境界条件
CHAR=AFFE_CHAR_MECA(MODELE=MODE,
                    DDL_IMPO=(_F(GROUP_MA=('kotei1','kotei2',), #2つの片持ち梁の固定端のグループ名
                                 DX=0,
                                 DY=0,
                                 DZ=0,),
                                 ),);

#(7)荷重条件(荷重が大きすぎると収束しない)
chr_no=AFFE_CHAR_MECA(MODELE=MODE,
                      FORCE_FACE=_F(GROUP_MA=('saika1','saika2',), #2つの片持ち梁の載荷面のグループ名
                                    #FZ=2.0,),);  # N/mm^2 : 引張の場合
                                    FY=0.03,),); # N/mm^2 : 曲げの場合

#以下、動的解析ではないが、他に適切な言葉が思い当たらないので、便宜上「時間」という言葉を使うので注意！

#(8)時間ごとの荷重の割合((時間,最大荷重の倍率)の並び)
depl_imp=DEFI_FONCTION(NOM_PARA='INST',VALE=(0,0,
                             0.5,1,              #time=0.5で最大荷重に到達して、
                             1,0,                #time=1.0でまた荷重0に戻る。
                             ),);

#(9)荷重増分
pas=DEFI_LIST_REEL(DEBUT=0.0,                    #time=0.0が開始
                   INTERVALLE=_F(JUSQU_A=1,      #time=1.0が終了
                                 NOMBRE=100,),); #time=0~1を100分割して計算する

#(10)収束計算の設定(海老さんが詳しい？)
pas2=DEFI_LIST_INST(DEFI_LIST=_F(METHODE='AUTO',
                                 LIST_INST=pas, #(9)の名前
                                 PAS_MINI=1e-05,
                                 PAS_MAXI=0.5,
                                 NB_PAS_MAXI=1000000,),
                    ECHEC=_F(EVENEMENT='ERREUR',
                             SUBD_PAS=4,
                             SUBD_PAS_MINI=1e-05, #許容誤差？
                             SUBD_METHODE='UNIFORME',),
                    ADAPTATION=_F(EVENEMENT='SEUIL',
                                  NB_INCR_SEUIL=2,
                                  NOM_PARA='NB_ITER_NEWTON',
                                  CRIT_COMP='LE',
                                  VALE_I=5,
                                  PCENT_AUGM=100,),);

#(11)静的な非線形解析の設定
RESU=STAT_NON_LINE(MODELE=MODE,
                   CHAM_MATER=MATE,
                   EXCIT=(_F(CHARGE=CHAR,),
                          _F(CHARGE=chr_no,
                             FONC_MULT=depl_imp,),),       #(8)の名前
                   COMP_INCR=_F(RELATION='VMIS_ISOT_TRAC', #弾塑性
                                DEFORMATION='PETIT',       #ひずみ・変位・回転が小さい場合
                                #DEFORMATION='GROT_GDEP',   #ひずみは小さいが、変位・回転が大きい場合
                                #DEFORMATION='PETIT_REAC',  #回転は小さいが、ひずみが大きい場合
                                #DEFORMATION='SIMO_MIEHE',  #ひずみ・変位・回転が大きい場合
                                #DEFORMATION='GDEF_LOG',    #ひずみ・変位・回転が大きい場合
                                GROUP_MA='bou1',),         #要素グループ名1
                   COMP_ELAS=_F(RELATION='ELAS_VMIS_TRAC', #弾性
                                DEFORMATION='PETIT',       #ひずみ・変位・回転が小さい場合
                                #DEFORMATION='GROT_GDEP',   #ひずみは小さいが、変位・回転が大きい場合
                                GROUP_MA='bou2',),         #要素グループ名2
                   INCREMENT=_F(LIST_INST=pas2,),          #(10)の名前
                   NEWTON=_F(REAC_ITER=1,),
                   CONVERGENCE=_F(RESI_GLOB_RELA=1e-6,
                                  ITER_GLOB_MAXI=500,),    #収束計算の最大反復数
                   ARCHIVAGE=_F(CHAM_EXCLU='VARI_ELGA',),);

#(12)要素解の出力
RESU=CALC_ELEM(reuse =RESU, #(11)の名前
               MODELE=MODE,
               CHAM_MATER=MATE,
               RESULTAT=RESU,
               OPTION=('SIEF_ELNO','SIEQ_ELNO','EPSI_ELNO','EPEQ_ELNO','EPSG_ELNO',),); #余計なものも色々ある

#(13)節点解の出力
RESU=CALC_NO(reuse =RESU, #(11)の名前
             RESULTAT=RESU, #(12)の名前
             OPTION=('SIEQ_NOEU','EPSI_NOEU','SIEF_NOEU','EPEQ_NOEU','FORC_NODA','REAC_NODA','EPSG_NOEU',), #余計なものも色々ある
             MODELE=MODE,
             CHAM_MATER=MATE,);

#(14)任意のグループの結果を.resuファイルに出力
#reacF=POST_RELEVE_T(ACTION=(_F(OPERATION='EXTRACTION',
#                               INTITULE='Reaction_Force_1',
#                               RESULTAT=RESU,
#                               NOM_CHAM='REAC_NODA',
#                               GROUP_NO='saika1',
#                               RESULTANTE='DZ',),
#                            _F(OPERATION='EXTRACTION',
#                               INTITULE='Reaction_Force_2',
#                               RESULTAT=RESU,
#                               NOM_CHAM='REAC_NODA',
#                               GROUP_NO='saika2',
#                               RESULTANTE='DZ',),),);
#
#IMPR_TABLE(TABLE=reacF,); #(14)の名前

#(15)出力の形式
IMPR_RESU(FORMAT='MED',
          UNITE=80,
          RESU=_F(RESULTAT=RESU,),); #(12),(13)の名前

FIN();
