有限元计算结构力学fortran程序
有限元计算结构力学fortran程序 计算结构力学程序 计算结构力学编程大作业 时间, 2007年6月 !!!**************************************************************************** !!! 关于程序的说明 !!!**************************************************************************** !一、功能: ! 1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平 ! 面桁架、梁、刚架及其组合结构的节点位移与杆端力; ! 2、可同时计算多种工况下的节点位移与杆端力。 !***************************************************************************** !****************************************************************************** ! ! 二、变量说明: ! NE——单元数; ! N——结构中自由度数; ! NJ——节点数; ! NS——特殊节点数，包括支座节点、主从节点(1节点不做主节点)、连接桁架的铰节点(没有转角); ! NAI——结构的单元截面类型数; ! MT——单元截面类型号; ! NL——荷载工况数; ! H——截面高度; ! E——弹性模量; ! JC——单元定位向量数组; ! X(NJ),Y(NJ)——节点的X，Y坐标值; ! JE(NE,2)——单元两端节点码数组; ! AI(NAI,2)——按单元类型顺序存放A与I，AI(I，1)—第I类单元的截面积，AI(I，2)—第I类单元的 ! 惯性矩; ! MT(NE)——单元所属单元类型号; ! JS(NS,4)——特殊节点信息，JS(I,1)—结点码;JS(I,2)，JS(I,3)，JS(I,4)—U，V，CETA约束信息， ! 有约束为1，没有约束为0;从节点某位移同主节点时位移时，该位移约束信息填主节点码; ! ! PJ(NP,3)——节点荷载信息数组;PJ(I，1)—节点力所在节点号;PJ(I，2)—节点力作用坐标方向: ! 坐标方向U，V，M分别为1，2，3; PJ(I，3)—节点力的大小(含正负号);U，V方向集中力时， ! 与坐标轴正向同向为正，M按右手法则为正;本程序推导过程取y轴向下为正。 ! ! PF(NF,4)——非节点荷载数组，并给出以下类型说明: ! 前6类型数据输法(梯形等可以用叠加法计算): ! PF(I,1)-单元码;PF(I,2)-类型;PF(I,3)-荷载大小;PF(I,4)-c值; ! 1——垂直于单元的均布力,大小为q，以坐标轴正向为正,c为荷载末端距i节点距离; ! 2——非节点集中力P，c为荷载距i节点距离; 1
计算结构力学程序 ! 3——非节点集中力距M，c为荷载距i节点距离，右手法则判正负; ! 4——三角形荷载，c为荷载距i节点距离，i端为0，距离i端c时力为q; ! j端为0的三角形，可按叠加法处理。 ! 5——沿杆轴向均布力,大小为q,c为荷载末端距i节点距离; ! 6——沿杆轴向集中力,大小为q,c为荷载末端距i节点距离; ! ! 从第7到第9类型(支座沉降)数据输法:PF(I,1)-单元码;PF(I,2)-类型;PF(I,3)-位移大小(含正负)，坐 标轴正向位为正，转角按右手法则;PF(I,4)-沉降所在的单元位移分量，i端为1-3，j端为4-6; ! ! 7——沿轴向支座沉降; ! 8——垂直于轴向支座沉降; ! 9——支座转动; ! 10——制造误差,PF(I,1)—制造误差所在单元，PF(I,2)-类型;PF(I,3)-误差大小(含正负),正负取决于消除 ! 误差时端点的运动方向，PF(I,4)—误差所在坐标号; ! 11——温度荷载，PF(I,1)—荷载所在单元，数据形式为:ElementNo.1 ,如2单元上有温度荷载，则PF(I,1)=2.1; ! PF(I,2)—温度变化值t1, PF(I,3)—温度变化值t2, PF(I,4)—材料线膨胀系数; ! ! TK(NN)——采用一维存储结构刚度矩阵，上半带元素(每列第一个非零元素到对角元); ! KD——主元地址数组，表示结构刚度矩阵的主元在TK中的序号，KD中最后一个数是TK中元素的总个数; ! JI——结构刚度矩阵上半带的非对角元素在TK中的地址，JI=KD(J)-J+I; ! JN(NJ,3)——结点位移分量编号数组，用于存放结点三个位移的位移分量号码， ! JN(I,1)，JN(I,2)，JN(I,2)-分别为结点I的U，V，CETA分量的位移分量(坐标)号码; ! ! P(N)——节点荷载列阵;在回代求位移时存放位移量; ! F(N)——求得的杆端力列阵; ! FO(6)——等效节点荷载列阵; ! ! !!!**************************************************************************************** !!!********************** 平面结构分析源程序内容 ************************************** !!!**************************************************************************************** PROGRAM PFF DIMENSION X(50),Y(50),JE(50,2),MT(30),AI(10,2),JS(20,4),PJ(50,3),PF(50,4),JN(50,3), & KD(150),TK(1000),P(150),F(6),H(50) DOUBLE PRECISION TK,P,F CHARACTER *200 TL OPEN(1,FILE='INDAT.DAT',STATUS='OLD') OPEN(2,FILE='OUTDAT.DAT',STATUS='NEW') READ(1,70) TL READ(1,70) TL READ(1,*)NE,NJ,NS,NAI,NL,E 2
计算结构力学程序 WRITE(2,10)NE,NJ,NS,NAI,NL,E 10 FORMAT(5X,'PLANE FRAME STRUCTURE ANALYSIS'/5X,'**********'//2X,'CONTROL PARAMETERS &OF STRUCTURE'/5X,'NE=',I2,8X,'NJ=',I2,8X,'NS=',I2,8X,'NAI=',I2,/5X,'NL=',I2,8X,'E=',E12.4) CALL INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H) !读入数据文件 CALL DJN(NJ,NS,JS,JN,N) !计算结构自由度数N，形成结点位移分量数组JN CALL ADE(NE,NJ,N,JE,JN,KD,NN) !形成主元地址数组KD(N) CALL SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK) !形成总刚，一维存储数组TK(NN) CALL UTDU3(TK,NN,KD,N) !对总刚进行UTDU分解，以用于解方程组 DO 20 LC=1,NL !对工况循环 READ(1,70)TL READ(1,70)TL READ(1,70)TL READ(1,*)NP,NF !读入工况信息 WRITE(2,30)LC,NP,NF 30 FORMAT(/2X,'LOAD DATA'/10X,'LOAD CASE=',I3/10X,'NP=',I3,8X,'NF=',I3) CALL NLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵P(N) CALL BACK3(TK,NN,P,N,KD,JN,NJ) !解方程组并输出结点位移，存放在数组P(N)中 WRITE(2,40) 40 FORMAT(//4X,'MEMBER-END FORCES OF ELEMENTS'/4X,'ELEMENT',13X,'N',17X,'V',17X,'M') DO 60 M=1,NE CALL MQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H) !计算单元杆端力，存放在数组F(6)中 WRITE(2,50)M,(F(I),I=1,6) !输出杆端力 50 FORMAT(/1X,I10,3X,'N1=',D12.4,3X,'V1=',D12.4,3X,'M1=',D12.4/14X,'N2=', &D12.4,3X,'V2=',D12.4,3X,'M2=',D12.4) 60 CONTINUE 20 CONTINUE 70 FORMAT(A) CLOSE(1) CLOSE(2) END SUBROUTINE INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H) !读入数据文件 DIMENSION X(NJ),Y(NJ),JE(NE,2),MT(NE),AI(NAI,2),JS(NS,4),H(NE) INTEGER NO READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)(NO,X(I),Y(I),I=1,NJ) READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)(NO,JE(I,1),JE(I,2),MT(I),H(I),I=1,NE) READ(1,70) TL
八节点平面等参元Fortran源程序
八节点平面等参元Fortran源程序C 这是一个采用平面四边形八节点等参单元的平面有限元分析程序。
PROGRAM PLANEFEMIMPLICIT REAL*8(A-H,O-Z)CHARACTER*80 LINECHAR,NEWLINECHARCOMMON A(30000),L(4000)COMMON /SOL/NPOIN,NELEM,NTYPE,NMATSOPEN(5,FILE='FEMDATA',STATUS='OLD')OPEN(6,FILE='FEMOUT',STATUS='UNKNOWN')C 以下进行的是从数据文件FEMDATA中读入数据文件主标题信息。
READ(5,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'输入')IF(LOCATECHAR.NE.0) THENLINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'ENDIFWRITE(6,5100) LINECHAR5000 FORMAT(A)5100 FORMAT(80('*')/A/80('*')/)C 以下进行的是先读入一行字符，然后从字符中读入网格单元总数NELEM。
READ(5,5000) LINECHARWRITE(6,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'=')LOCATECHAR=LOCATECHAR+1NEWLINECHAR=LINECHAR(LOCATECHAR:80)READ(NEWLINECHAR,5200) NELEM5200 FORMAT(I5)C 以下进行的是先读入一行字符，然后从字符中读入单元节点总数NPOIN。
READ(5,5000) LINECHARWRITE(6,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'=')LOCATECHAR=LOCATECHAR+1NEWLINECHAR=LINECHAR(LOCATECHAR:80)READ(NEWLINECHAR,5200) NPOINC 以下进行的是先读入一行字符，然后从字符中读入问题类型编号NTYPE。
Fortran程序设计(第六章-循环结构(下))
【5】为了计算并输出n!,其中n从键盘输入，下 列各FORTRAN程序中正确的是： D) C) READ(*,*)N READ(*,*)N K=1 K=1 S=1.0 S=1.0 10 IF(K.LE.N) THEN 10 IF(K.LE.N)THEN K=K+1 S=S*K S=S*K K=K+1 GOTO 10 GOTO 10 END IF END IF WRITE(*,*)’S=’,S WRITE(*,*)’S=’,S END END
练习：
1、下面关于DO循环的规定，错误的是______。 (A) DO循环的循环控制变量不能在循环体内赋值 (B) DO循环的控制变量表达式,终值和步长可以是整型和实型 (C) DO循环是当型循环 (D) DO循环的循环控制变量不能是双精度型 2、 有如下循环入口语句 INTEGER：：I DO I=-0.5,-0.5,-1.0 该循环的执行次数为____________。 (A) 0 (B) 1 (C) 出错 (D) 无限
Do结构可以有多重嵌套，这里介绍二重嵌套的执行过程。 对于多重嵌套，其基本原理相同。
1．当控制进入到外层DO结构后，先计算出外层DO结构的 循环次数Ri，外层循环变量得到初值。 2．若Ri<0，则结束外循环的执行，当然也不能进入内循环； 若Ri>0，执行外层结构的DO块内的语句。 3．当遇到内层DO语句时，控制进入内层DO结构；先算出 内层DO结构的循环次数Rj，内层循环变量得到初值。
第六章 循环结构（下）
§6.5 DO结构嵌套
§6.6 隐含DO循环 §6.7 程序举例
6.5 DO结构嵌套
一个DO结构循环体内可以包含另一个DO循环结构，
这就是DO结构循环嵌套。 注意： 1 内循环必须完全嵌套在外循环体内，不能相互交叉。
Fortran程序设计第11章 基本计算(三)循环控制结构
第11章 基本计算(三)循环控制结构 上章讨论的控制结构的特点是通过判别条件来对结构内的块进行选择，所对应的算法结构，最简单的例子，就是解一元二次方程，在输入方程所以的参数值之后，需要首先计算一个判别式，然后根据判别式的值，再选择使用哪个公式，也就是计算的途径，才能够给出最终的解。
在本章所讨论的控制结构的特点则是针对结构内的块进行多次的重复运算，每完成一次运算，都判别一下是否需要把此次运算结果作为输入，再进行一次运算。
这种控制结构所对应的算法结构，一个最简单的例子，就是求级数的部分和。
我们知道对于具有通项表达式的级数，求它的部分和的每一项，总是需要进行同样的运算过程，如果使用按照序列形式排列的程序结构，那么需要计算多少项，就需要写下多少条语句，把它们顺序排列下来，才能做到程序走一遍即完成计算。
这样的算法设计显然是没有利用运算过程里所表现的循环结构，如果使用一种控制结构与循环过程对应，让程序的运行能够重复表示通项公式的表达式计算，就能够用一个表达式代替所有项的表达式，显然更加合理。
【例11-1】 设有一个级数：311/ni x i ==∑如果要求级数在N=K 时的值，如果一定要使用序列结构的程序，那么在程序当中肯定会出现如下K 个表达式顺序排列的情形：…I=1SUM=1/I**3I=2SUM=SUM+1/I**3I=3SUM=SUM+1/I**3…I=KSUM=SUM+1/I**3由于这K 个表达式是一样的，因此如果使用如下的一个控制结构，只需要使用一个表达式赋值语句，就可以表示整个循环运算过程：…SUM=0.0DO I=1，KSUM= SUM+1/I**3END DO上面的控制结构，就是本章所要讨论的DO结构。
实际上这是一种基本的计算过程，为很多重要算法的实现提供了基础。
FORTRAN语言提供用来进行循环控制的主要结构就是DO结构，因此本章主要讨论的就是DO结构。
最后还会简略的讨论有关分支转移的实现问题，尽管它们不属于循环控制，但是由于在现代结构性编程风格的要求下，这种分支转移是一种过时的方式，因此简略地附加在本章后面。
Matlab 有限元法计算分析程序编写
6) M函数文件 与命令文件不同，函数文件从外界只能看到传给它的输入 量和送出来的计算结果，而内部运作是看不见的。它的特 点是 (1)从形式上看，与命令文件不同，函数文件的第一行总是以 “function”引导的“函数申明行”。 (2)从运行上看，与命令文件运行不同，每当函数文件运行， MATLAB就会专门为它开辟临时工作空间，所有中间变量 都存放在函数工作空间中，当执行完文件最后一条指令和 遇到return时，就结束该函数文件的执行，同时该临时函数 工作空间及其所有的中间变量立即被清除。 (3)对于函数文件中的变量，如果不作特别说明，默认为临时 局部变量，这些临时变量就存放在函数的临时工作空间中， 当函数结束时他们被立即清除。与之相对应的是全局变量， 他们是通过global指令进行特别申明，这些全局变量可被几 个不同的函数共享。 • 函数文件的编辑也可用MATLAB editor/debugger。
有限元法计算分析程序编写
结构参数输入，包括
1）节点坐标值 2）单元类型以及连接信息 3）各单元的弹性模量、截面积（厚度）等 4）荷载形式以及作用位置、作用方向、荷载值 5）约束条件 6）输出信息
m j
对节点和单元分别编号 每个节点的自由度根据 节点号计算得到
i
y
o
x
计算结构的刚度矩阵
对各单元作如下的计算 a)计算单元刚度矩阵 b)计算坐标转换矩阵（如果需要） c)作坐标转换计算（如果需要） d)按自由度顺序叠加到总刚度矩阵中
MATLAB的使用方法
1) 最简单的计算器使用法 求[12＋2×(7－4)]÷32的算术运算结果 （1）用键盘在MATLAB指令窗中输入一下内容 (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按【Enter】键，该指令被执行 （3）在指令执行后，MATLAB指令窗中将显示一下内容 ans = 2 [说明] 加 + 减 乘 * 除 / 或 \ (这两个符号对于数组有不同的含义) 幂 ^ “ans”是answer的缩写，其含义是运算答案，它是MATLAB的一个默 认变量
空间框架静力计算源程序(程序结构力学大作业)
主要为一些工具函数 变量定义，单元属性分析 矩阵求解 整体控制模块
1：Lxz_Tools.f90 module Lxz_Tools implicit none integer (kind(1)),parameter ::ikind=(kind(1)) integer (kind(1)),parameter ::rkind=(kind(0.D0)) real (rkind), parameter :: Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0, & & Four=4.D0,Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0,Nine=9.D0, & & Ten=10.D0 contains function matinv(A) result (B) real(rkind) ,intent (in)::A(:,:) !real(rkind) , allocatable::B(:,:) real(rkind) , pointer::B(:,:) integer(ikind):: N,I,J,K real(rkind)::D,T real(rkind), allocatable::IS(:),JS(:) N=size(A,dim=2) allocate(B(N,N)) allocate(IS(N));allocate(JS(N)) B=A do K=1,N D=0.0D0 do I=K,N do J=K,N if(abs(B(I,J))>D) then D=abs(B(I,J)) IS(K)=I JS(K)=J
平面三角形3结点有限元程序
一、平面三角形３结点有限元程序1、程序名：，FEM3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 5.0编制而成，输入数据全部采用自由格式。
3、程序流程及框图图１-1 程序流程图图１-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移 ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。
4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。
数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR ——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束力的结点数目。
有限元计算
有限元计算
有限元计算（Finite Element Analysis，FEA）是一种数值计算
方法，用于计算复杂结构在外部载荷作用下的响应。
其基本思想是将结构分割成有限数量的元素，并在每个元素上进行力和位移的计算，然后将结果整合到整个结构中。
有限元计算可以应用于各种不同类型的工程领域，如航空航天、汽车工程、建筑工程、机械工程等。
它可以用来计算结构在极限载荷下的强度和稳定性，分析结构的动态响应，以及进行优化设计等。
有限元计算的基本步骤包括建立有限元模型、定义材料力学特性、定义边界条件和荷载、进行计算求解、分析结果并评估结构安全性。
在此过程中，需要使用特定的有限元软件和计算机资源来进行计算。
(整理)Fotran90版—平面刚架有限元分源程序代码.
Fotran90版—平面刚架有限元分源程序代码program mainreal,allocatable::ks(:,:)allocatable lnd(:,:)allocatable crd(:,:)allocatable ea(:)allocatable ei(:)allocatable jcs(:,:)allocatable pj(:,:)allocatable bl(:)allocatable p(:)open (5,file="inputdates.in")read (5,*) ne,nj,ns,npjnj3=3*njallocate (lnd(ne,2))allocate (crd(nj,2))allocate (ea(ne))allocate (ei(ne))allocate (jcs(ns,4))allocate (pj(npj,4))allocate (bl(ne))allocate (ks(nj3,nj3))allocate (p(nj3))write (*,"(1x,'plane fram structural analysis'//)")write (*,"(1x,'structural parameters'/)")write (*,"(/1x,'total number of')")write (*,"(1x,'element=',i5/1x,'joints=',i5/1x,'constructed joints=',i5/1x,'loads=',i5/)")&ne,nj,ns,npjcall readin(ne,nj,ns,npj,lnd,crd,ea,ei,jcs,pj,bl)call formf(nj3,npj,pj,p)call cks(nj3,ne,ea,ei,bl,lnd,crd,nj,ks)call dealbc(ns,jcs,nj3,ks,p)call solve(nj,nj3,ks,p)close (5)stopend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine readin(ne,nj,ns,npj,lnd,crd,ea,ei,jcs,pj,bl)dimensionlnd(ne,2),crd(nj,2),ea(ne),ei(ne),jcs(ns,4),pj(npj,4),bl(ne)read (5,*) ((lnd(i,j),j=1,2),i=1,ne)write (*,"(/1x,'element dates',/1x,'element',4x,'conection',8x)") write (*,"(1x,i5,2x,i5,3x,'to',i5,3x)") (i,(lnd(i,j),j=1,2),i=1,ne)read (5,*) ((crd(i,j),j=1,2),i=1,nj)write (*,"(1x,'nodelcoordinates'/3x,'node',6x,'x-coordinates',7x,'y-coordinates')")write (*,"(1x,i5,5x,f10.4,10x,f10.4)") (i,(crd(i,j),j=1,2),i=1,nj)read (5,*) (ea(i),i=1,ne)read (5,*) (ei(i),i=1,ne)write (*,"(/1x,'materialparameters',/1x,'element',11x,'ea',13x,'ei')")write (*,"(1x,i5,5x,2e15.6)") (i,ea(i),ei(i),i=1,ne)read (5,*) ((jcs(i,j),j=1,4),i=1,ns)write(*,"(/1x,'constrained nodes',/3x,'nodes',1x,'X',4x,'Y',4x,'R')") write (*,"(4i5)") ((jcs(i,j),j=1,4),i=1,ns)read (5,*) ((pj(i,j),j=1,4),i=1,npj)write (*,"(/1x,'joints oflonds',/1x,'joint',5x,'PX',8x,'PY',8x,'Mxy')")write (*,"(1x,f5.0,3f10.4)") ((pj(i,j),j=1,4),i=1,npj)do ie=1,nei=lnd(ie,1)j=lnd(ie,2)dx=crd(j,1)-crd(i,1)dy=crd(j,2)-crd(i,2)bl(ie)=sqrt(dx**2+dy**2)end dowrite (*,"(/1x,'the length of the elements',/1x,'elementnumers',5x,'numembers length')")write (*,"(1x,i5,10x,f10.4)") (i,bl(i),i=1,ne)returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine turn(ne,lnd,crd,bl,t,ie,nj) dimension lnd(ne,2),crd(nj,2),bl(ne),t(6,6)i1=lnd(ie,1)j1=lnd(ie,2)dx=crd(j1,1)-crd(i1,1)dy=crd(j1,2)-crd(i1,2)si=dy/bl(ie)co=dx/bl(ie)do i=1,6t(i,1:6)=0.0end dot(1,1)=cot(1,2)=sit(2,1)=-sit(2,2)=cot(3,3)=1.0do i=1,3do j=1,3t(i+3,j+3)=t(i,j)end doend do!write (*,"(/1x,'the dates of t')")!write (*,"(/1x,6f10.4)") ((t(i,j),j=1,6),i=1,6) returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine stif(ne,ea,ei,bl,kd,ie)dimension ea(ne),ei(ne),bl(ne)real kd(6,6)a1=ea(ie)e1=ei(ie)s=bl(ie)do i=1,6kd(i,1:6)=0.0end dokd(1,1)=a1/skd(2,2)=12.0*e1/s**3kd(3,2)=-6.0*e1/s**2kd(3,3)=4.0*e1/skd(4,1)=-kd(1,1)kd(4,4)=kd(1,1)kd(5,2)=-kd(2,2)kd(5,3)=-kd(3,2)kd(6,3)=2.0*e1/skd(6,5)=-kd(3,2)kd(6,6)=kd(3,3)do i=1,6do j=1,ikd(j,i)=kd(i,j)end doend doreturnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine stie(ne,ea,ei,bl,kd,lnd,crd,t,ke,ie,nj)dimension ea(ne),ei(ne),bl(ne),lnd(ne,2),crd(nj,2),t(6,6),ek(6,6) real kd(6,6),ke(6,6)call stif(ne,ea,ei,bl,kd,ie)call turn(ne,lnd,crd,bl,t,ie,nj)do i=1,6do j=1,6ek(i,j)=0.0do k=1,6ek(i,j)=ek(i,j)+kd(i,k)*t(k,j)end doend doend dodo i=1,6do j=1,6ke(i,j)=0.0do k=1,6ke(i,j)=ke(i,j)+t(k,i)*ek(k,j)end doend doend doreturnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine formf(nj3,npj,pj,p)dimension p(nj3),pj(npj,4)p=0.0j=pj(i,1)p(3*j-2)=p(3*j-2)+pj(i,2)p(3*j-1)=p(3*j-1)+pj(i,3)p(3*j)=p(3*j)+pj(i,4)end do!write(*,"(/1x,f10.4)") (p(i),i=1,nj3)returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine cks(nj3,ne,ea,ei,bl,lnd,crd,nj,ks) dimension ea(ne),ei(ne),bl(ne),lnd(ne,2),crd(nj,2),& t(6,6),ek(6,6)real ks(nj3,nj3),kd(6,6),ke(6,6)ks(1:nj3,1:nj3)=0.0do ie=1,necall stie(ne,ea,ei,bl,kd,lnd,crd,t,ke,ie,nj)do i=1,2do ii=1,3ir=3*(i-1)+iiiw=3*(lnd(ie,i)-1)+iido j=1,2do jj=1,3jr=3*(j-1)+jjjw=3*(lnd(ie,j)-1)+jjks(iw,jw)=ks(iw,jw)+ke(ir,jr)end doend doend doend doend do!write (*,"(/1x,'dates of ks')")!write (*,"(/1x,9e10.3)") ((ks(i,j),j=1,nj3),i=1,nj3) returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine dealbc(ns,jcs,nj3,ks,p)dimension jcs(ns,4),p(nj3)real ks(nj3,nj3)i1=jcs(i,1)do j=1,3j1=jcs(i,j+1)if (j1==0)thenelseiw=3*(i1-1)+jdo k=1,nj3if (k==iw)thenks(iw,iw)=1.0elseks(iw,k)=0.0ks(k,iw)=0.0end ifend dop(iw)=0.0end ifend doend do!write (*,"(/1x,'the dates of bc')")!write (*,"(/1x,9e15.6,5x,e15.6)") ((ks(i,j),j=1,nj3),p(i),i=1,nj3) returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine solve(nj,nj3,ks,p)dimension p(nj3)real ks(nj3,nj3)do k=1,nj3-1do i=k+1,nj3c=ks(i,k)/ks(k,k)do j=k,nj3ks(i,j)=ks(i,j)-c*ks(k,j)end dop(i)=p(i)-c*p(k)end doend dop(nj3)=p(nj3)/ks(nj3,nj3)do i=nj3-1,1,-1do j=i+1,nj3p(i)=p(i)-ks(i,j)*p(j)end dop(i)=p(i)/ks(i,i)end dowrite (*,"(/1x,'displacement of the joints',/1x,'jointsnumbers',3x,'X-dis',10x,'Y-dis',10x,'Z-rot')")write (*,"(1x,i5,5x,3e15.6)") (i,p(3*i-2),p(3*i-1),p(3*i),i=1,nj) returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!结果输出如下：。
fortran常用算法程序集
Fortran常用算法程序集简介Fortran是一种面向科学和工程计算的编程语言，通常用于数值计算和数据分析。
它有着强大的数学计算能力和高性能，被广泛应用于科学计算、工程仿真、天气预报等领域。
本文将介绍一些常用的Fortran算法程序集，包括数值积分、矩阵运算、排序算法等。
数值积分数值积分是求解定积分的一种方法，用于计算曲线下面积、求解微分方程等。
Fortran提供了一些常用的数值积分算法，如梯形法则、辛普森法则等。
梯形法则梯形法则是数值积分中最简单的算法之一，基本思想是将曲线下面积近似为一系列梯形的和。
下面是使用Fortran编写的梯形法则算法示例：! 梯形法则real function trapezoidal_rule(f, a, b, n)real, external:: freal:: a, b, n, hreal:: x, suminteger:: ih = (b - a) / nsum= f(a) + f(b)do i =1, n-1x = a + i * hsum=sum+2* f(x)end dotrapezoidal_rule = (h/2) *sumend function trapezoidal_rule辛普森法则辛普森法则是一种更精确的数值积分算法，基于多项式插值的思想。
它将曲线分成若干小段，每段近似为一个二次函数，然后对每个二次函数进行积分。
下面是使用Fortran编写的辛普森法则算法示例：! 辛普森法则real function simpsons_rule(f, a, b, n)real, external:: freal:: a, b, n, hreal:: x, sum1, sum2integer:: ih = (b - a) / nsum1 = f(a) + f(b)sum2 =0do i =1, n-1, 2x = a + i * hsum2 = sum2 +4* f(x)end dodo i =2, n-2, 2x = a + i * hsum2 = sum2 +2* f(x)end dosimpsons_rule = (h/3) * (sum1 + sum2)end function simpsons_rule矩阵运算矩阵运算是科学计算中常用的一个重要环节，Fortran提供了丰富的矩阵运算库，包括矩阵乘法、矩阵转置、矩阵求逆等。
有限元基础及程序详解
d e dr
That is, the component of in the direction of e gives the rate of change of
in that direction (the directional derivative). In particular, the components of
in the coordinate directions e i are given by
d ei dr in the ei direction xi
Therefore, the Cartesian components of are
, that is, xi
d (r dr) - (r) dr
If dr denote the magnitude of dr , and e the unit vector in the direct ion of
dr ( Note: e = dr dr ). Then the above equation gives, for dr in the e direction.
div v tr(v)
In Cartesian coordinates, this gives
div v
v1 v2 v3 vm x1 x2 x3 xm
Let T(r ) be a tensor field. The divergence of T(r ) is defined to be a vector field, denoted by divT , such that for any vector a
T( a + b) = Ta + Tb