fortran 悬臂梁模型代码

### Fortran 编写悬臂梁模型

以下是使用 Fortran 实现的一个简单悬臂梁有限元分析的示例代码。该代码主要用于求解线弹性条件下悬臂梁的变形情况：

program cantilever_beam_fem
    implicit none
    
    integer, parameter :: ndof = 2      ! 每个节点自由度数
    integer, parameter :: nnode = 5     ! 节点总数
    integer, parameter :: nelem = 4     ! 单元总数
    real*8, dimension(ndof*nnode) :: d  ! 位移向量
    real*8, dimension(nnode, 2) :: coord! 坐标数组
    real*8, dimension(nelem, 2) :: conn ! 连接关系表
    real*8, dimension(ndof*nnode, ndof*nnode) :: k_global ! 全局刚度矩阵
    real*8, dimension(ndof*nnode) :: f_external ! 外力载荷向量
    real*8 E, A, L, I                    ! 材料属性和几何尺寸变量
    integer i, j                          ! 循环计数器
    
    ! 初始化材料参数和截面特性 (单位制需统一)
    E = 70e9          ! 弹性模量 Pa
    A = 1.e-4         ! 截面积 m^2
    I = 1./12 * A**3 / 1.           ! 抗弯惯性矩 m^4
    L = 1.            ! 总长度 m
    
    ! 定义网格划分
    do i=1,nnode
        coord(i, :) = [(i-1)*L/(nnode-1), 0.d0]
    enddo
    
    ! 设置单元连接关系
    do i=1,nelem
        conn(i,:) = [i,i+1]
    enddo
    
    call assemble_stiffness_matrix(k_global, E, I, L/nelem, conn, nnode, nelem)

    ! 施加边界条件并组装外力向量
    forall(j=1:nnode*ndof) 
        if(mod(j,2)==1 .and. j<=ndof*(nnode/2)) then
            f_external(j)=0.
        elseif((j>=(nnode-1)*ndof).and.(mod(j,2)/=1))then
            f_external(j)=-1000.*E*A*L/nnode   ! 应用端部集中力
        endif
    endforall
    
    ! 解方程组 Kd=F 并输出结果
    call solve_linear_system(k_global,f_external,d)
    
contains

subroutine assemble_stiffness_matrix(K,E,I,L,e_conn,n_node,n_elem)
    use iso_fortran_env
    implicit none
    real(kind=real64), intent(out)::K(:,:)
    real(kind=real64),intent(in)::E,I,L
    integer,intent(in):: e_conn(:,:),n_node,n_elem
    real(kind=real64),dimension(4,4)::k_local
    integer:: elem,node,mapping(4)
    
    ! 计算局部刚度矩阵
    k_local = reshape([&
        & 12*I/L**3 , 6*I/L**2 ,-12*I/L**3 , 6*I/L**2,&
        & 6*I/L**2  , 4*I/L   ,-6*I/L**2  , 2*I/L  ,&
        &-12*I/L**3,-6*I/L**2 , 12*I/L**3 ,-6*I/L**2,&
        & 6*I/L**2  , 2*I/L   ,-6*I/L**2  , 4*I/L], &
        &[4,4])
        
    ! 组装全局刚度矩阵
    K = 0._real64
    do elem=1,n_elem
        mapping=[(conn(elem,j)-1)*2+j,j=1,2]*2-1
        do node=1,4
            do i=node,4
                K(mapping(node),mapping(i))=K(mapping(node),mapping(i))+k_local(node,i)
                if(node/=i) K(mapping(i),mapping(node))=K(mapping(node),mapping(i))
            enddo
        enddo
    enddo
    
end subroutine assemble_stiffness_matrix


subroutine solve_linear_system(A,b,x)
    use lapack_driver
    implicit none
    real(kind=real64),intent(inout)::A(:,:),b(:)
    real(kind=real64),allocatable::work(:)
    integer::ipiv(size(b)),info,lwork,size_b
    size_b=size(b)
    allocate(work(max(1,3*size_b)))
    lwork=max(1,3*size_b)
    call DGETRF(size_b,size_b,A,size_b,ipiv,info)
    if(info==0)call DGETRS('N',size_b,1,A,size_b,ipiv,b,size_b,info)
    deallocate(work)
end subroutine solve_linear_system

end program cantilever_beam_fem

此代码实现了以下几个功能：
- 创建了一个简单的二维平面内的离散化悬臂梁；
- 使用四阶多项式插值来构建每个梁单元的局部刚度矩阵[^1]；
- 将所有单元贡献累加以形成整体结构系统的总刚度矩阵；
- 对于固定端施加约束，在自由端加载垂直方向上的集中力。