[讨论]

纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之c

2022-11-25 10:46:03 683 控制程序

0 !*********************************************************** !********************************************************* ! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之c !********************************************************* ! SINOMACH ! ! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之c ! ! 周勇 20221125 !********************************************************** !c $$$ Run, Output: !********************************************************** !********************************************************** Subroutine p_write(p_file_name, node_max, node_num, neqn, c, indx, node_x, node_y) !******************************************************************72 ! !c P_WRITE writes the pressures values for a given timestep to a file. ! Implicit None Integer node_max Integer node_num Integer neqn Double Precision c(neqn) Character (Len=) p_file_name Integer i Integer indx(node_max, 3) Integer node Double Precision node_x(node_num) Double Precision node_y(node_num) Double Precision p Open (Unit=1, File=p_file_name, Err=10) Do node = 1, node_num `i = indx(node, 3) If (0==i) Then Else If (0<i) Then p = c(i) Write (1, '(g14.6)') p Else If (i<0) Then p = 0.0D+00 Write (1, '(g14.6)') p End If` End Do Close (Unit=1) Write (, '(a)') ' ' Write (, '(a)') 'P_WRITE:' Write (, '(a)') ' Wrote the pressure output file:' Write (, '(a)') ' ' Write (, '(4x,a)') p_file_name Return 10 Continue Write (, '(a)') ' ' Write (, '(a)') 'P_WRITE - Fatal errorc' Write (, '(a)') ' Could not open the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') p_file_name Stop End Subroutine p_write Subroutine quad_a_set(maxel, nelemn, nquad, area, xm, ym) !*******************************************************************72 ! !c QUAD_A_SET sets quadrature information for the assembly routine. ! Implicit None Integer maxel Integer nquad Double Precision area(maxel) Integer element Integer nelemn Double Precision xm(maxel, nquad) Double Precision ym(maxel, nquad) Do element = 1, nelemn xm(element, 1) = 0.5D0 xm(element, 2) = 0.5D0 xm(element, 3) = 0.0D0 End Do Do element = 1, nelemn ym(element, 1) = 0.0D0 ym(element, 2) = 0.5D0 ym(element, 3) = 0.5D0 End Do Do element = 1, nelemn area(element) = 0.5D0 End Do Return End Subroutine quad_a_set Function r8_huge() !*****************************************************************72 ! !c R8_HUGE returns a "huge" double precision number. ! Implicit None Double Precision r8_huge r8_huge = 1.0D+30 Return End Function r8_huge Function refbsp(x, y, iq) !******************************************************************72 ! !c REFBSP evaluates a linear basis functions on the reference triangle. ! Implicit None Integer iq Double Precision refbsp Double Precision x Double Precision y If (iq==1) Then refbsp = 1.0D+00 - x - y Else If (iq==2) Then refbsp = x Else If (iq==3) Then refbsp = y Else Write (, '(a)') ' ' Write (, '(a)') 'REFBSP - Fatal errorc' Write (, '(a)') ' Illegal input value of IQ = ', iq Stop End If Return End Function refbsp Subroutine refqbf(x, y, in, bb, bx, by) !********************************************************************72 ! !c REFQBF evaluates quadratic basis functions on the reference triangle. ! Implicit None Double Precision bb Double Precision bx Double Precision by Integer in Double Precision x Double Precision y If (in==1) Then bb = (1.D0-x-y)(1.D0-2.D0x-2.D0y) bx = -3.D0 + 4.D0x + 4.D0y by = -3.D0 + 4.D0x + 4.D0y Else If (in==2) Then bb = x(2.D0x-1.D0) bx = 4.D0x - 1.D0 by = 0.D0 Else If (in==3) Then bb = y(2.D0y-1.D0) bx = 0.D0 by = 4.D0y - 1.D0 Else If (in==4) Then bb = 4.D0x(1.D0-x-y) bx = 4.D0(1.D0-2.D0x-y) by = -4.D0x Else If (in==5) Then bb = 4.D0xy bx = 4.D0y by = 4.D0x Else If (in==6) Then bb = 4.D0y(1.D0-x-y) bx = -4.D0y by = 4.D0(1.D0-x-2.D0y) Else bb = 0.0D+00 bx = 0.0D+00 by = 0.0D+00 End If Return End Subroutine refqbf Subroutine setgrd(xc, yc, area, xm, ym, region_density_x,& region_density_y, region_x, region_y, node, indx, nlband,& nuband, nband, nelemn, np, nnodes, nuk, nquad, neqn, maxnd, maxel) !*******************************************************************72 ! !c SETGRD sets up the grid for the problem. ! Implicit None Integer maxel Integer maxnd Integer nnodes Integer nquad Integer nuk Double Precision area(maxel) Integer indx(maxnd, nuk) Integer nband Integer nelemn Integer neqn Integer nlband Integer node(maxel, nnodes) Integer np Integer nuband Integer region_density_x(3) Integer region_density_y(3) Double Precision region_x(4) Double Precision region_y(4) Double Precision xc(maxnd) Double Precision xlngth Double Precision xm(maxel, nquad) Double Precision yc(maxnd) Double Precision ylngth Double Precision ym(maxel, nquad) ! ! Determine NP, the number of nodes. ! Call hcell_node_count(region_density_x, region_density_y, np) ! ! Assign coordinates to the nodes, XC and YC. ! Call hcell_node_xy(region_density_x, region_density_y, np, region_x,& region_y, xc, yc) ! ! Determine NELEMN, the number of elements. ! Call hcell_element_count(region_density_x, region_density_y, nelemn) ! ! Assign nodes to elements in NODE. ! Call hcell_element_node(region_density_x, region_density_y, maxel,& nelemn, nnodes, node) ! ! Determine NEQN, the number of unknowns. ! Call hcell_dof_count(region_density_x, region_density_y, neqn) ! ! Assign degrees of freedom in INDX. ! Call hcell_dof_set(region_density_x, region_density_y, maxnd, np, indx) ! ! For the assembly routine, determine the quadrature data ! NQUAD, AREA, XM and YM. ! Call quad_a_set(maxel, nelemn, nquad, area, xm, ym) ! ! Compute the bandwidths NLBAND, NUBAND and NBAND. ! Call element_node_bandwidth(maxnd, maxel, nnodes, nelemn, node,& neqn, np, indx, nlband, nuband, nband) Return End Subroutine setgrd Subroutine timestamp() !******************************************************************72 ! !c TIMESTAMP prints out the current YMDHMS date as a timestamp. ! Implicit None Character (8) ampm Integer d Character (8) date Integer h Integer m Integer mm Character (9) month(12) Integer n Integer s Character (10) time Integer y Save month Data month/'January ', 'February ', 'March ', 'April ',& 'May ', 'June ', 'July ', 'August ', 'September',& 'October ', 'November ', 'December '/ Call date_and_time(date, time) Read (date, '(i4,i2,i2)') y, m, d Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm If (h<12) Then ampm = 'AM' Else If (h==12) Then If (n==0 .And. s==0) Then ampm = 'Noon' Else ampm = 'PM' End If Else h = h - 12 If (h<12) Then ampm = 'PM' Else If (h==12) Then If (n==0 .And. s==0) Then ampm = 'Midnight' Else ampm = 'AM' End If End If End If Write (, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,& month(m), y, h, ':', n, ':', s, '.', mm, ampm Return End Subroutine timestamp Subroutine trans(it, xq, yq, det, pj11, pj21, pj12, pj22, xc, yc, node, maxel) !********************************************************************72 ! !c TRANS transforms data between the reference and physical elements. ! Implicit Double Precision (A-H, O-Z) Double Precision xc() Double Precision yc() Integer node(maxel, ) i1 = node(it, 1) i2 = node(it, 2) i3 = node(it, 3) i4 = node(it, 4) i5 = node(it, 5) i6 = node(it, 6) x1 = xc(i1) y1 = yc(i1) x2 = xc(i2) y2 = yc(i2) x3 = xc(i3) y3 = yc(i3) x4 = xc(i4) y4 = yc(i4) x5 = xc(i5) y5 = yc(i5) x6 = xc(i6) y6 = yc(i6) ! ! Compute partial derivatives at point (xq,yq) ! f1x = x1(-3.D0+4.D0xq+4.D0yq) + x2(4.D0xq-1.D0) +& x44.D0(1.D0-2.D0xq-yq) + x54.D0yq + x64.D0(-yq) f1y = x1(-3.D0+4.D0xq+4.D0yq) + x3(4.D0yq-1.D0)& x44.D0(-xq) + x54.D0xq + x64.D0(1.D0-xq-2.D0yq) f2x = y1(-3.D0+4.D0xq+4.D0yq) + y2(4.D0xq-1.D0) & y44.D0(1.D0-2.D0xq-yq) + y54.D0yq + y64.D0(-yq) f2y = y1(-3.D0+4.D0xq+4.D0yq) + y3(4.D0yq-1.D0) & y44.D0(-xq) + y54.D0xq + y64.D0(1.D0-xq-2.D0yq) ! ! Compute determinant of transformation evaluated at point (xq,yq) ! det = f1xf2y - f1yf2x ! ! Compute j11, j22, j21, j22 ! pj11 = f2y/det pj12 = -f2x/det pj21 = -f1y/det pj22 = f1x/det det = dabs(det) Return End Subroutine trans Function ubdry(iuk, ip, xc, yc) !********************************************************************72 ! !c UBDRY evaluates the boundary conditions. ! Implicit None Integer ip Integer iuk Double Precision ubdry Double Precision xc() Double Precision yc() ubdry = 0.0D+00 If (xc(ip)==0.0D+00) Then If (iuk==1) Then ubdry = 16.0D+00(yc(ip)-0.5D+00)(1.0D+00-yc(ip)) End If End If Return End Function ubdry Subroutine uv_read(uv_file_name, neqn, uold) !*******************************************************************72 ! !c UV_READ returns an initial value for the solution coefficient vector. ! Implicit None Integer neqn Character (Len=) uv_file_name Integer i Double Precision uold(neqn) Open (Unit=1, File=uv_file_name, Err=10) Do i = 1, neqn Read (1, ) uold(i) End Do Close (Unit=1) Write (, '(a)') ' ' Write (, '(a)') 'UV_READ:' Write (, '(a)') ' Read the file containing' Write (, '(a)') ' the initial solution values:' Write (, '(a)') ' ' Write (, '(4x,a)') uv_file_name Return 10 Continue Write (, '(a)') ' ' Write (, '(a)') 'UV_READ - Fatal errorc' Write (, '(a)') ' Could not open the file containing' Write (, '(a)') ' the initial solution values:' Write (, '(a)') ' ' Write (, '(4x,a)') uv_file_name Stop End Subroutine uv_read Subroutine uv_write(uv_file_name, node_max, node_num, neqn, c, indx,& node_x, node_y) !*******************************************************************72 ! !c UV_WRITE writes the velocity values for a given timestep to a file. ! Implicit None Integer node_max Integer node_num Integer neqn Double Precision c(neqn) Character (Len=) uv_file_name Integer i Integer indx(node_max, 3) Integer node Double Precision node_x(node_num) Double Precision node_y(node_num) Double Precision u Double Precision v Open (Unit=1, File=uv_file_name, Err=10) Do node = 1, node_num `i = indx(node, 1) If (0<i) Then u = c(i) Else u = 0.0D+00 End If i = indx(node, 2) If (0<i) Then v = c(i) Else v = 0.0D+00 End If Write (1, '(2g14.6)') u, v` End Do Close (Unit=1) Write (, '(a)') ' ' Write (, '(a)') 'UV_WRITE:' Write (, '(a)') ' Wrote the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') uv_file_name Return 10 Continue Write (, '(a)') ' ' Write (, '(a)') 'UV_WRITE - Fatal errorc' Write (, '(a)') ' Could not open the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') uv_file_name Stop End Subroutine uv_write Subroutine xy3_write(file_name, maxnd, node_num, indx, node_x, node_y) !********************************************************************72 ! !c XY3_WRITE writes the coordinates of nodes of 3-node triangles to a file. ! Implicit None Integer maxnd Integer node_num Character (Len=) file_name Integer indx(maxnd, 3) Integer node Double Precision node_x(node_num) Double Precision node_y(node_num) Open (Unit=1, File=file_name, Err=10) Do node = 1, node_num If (0/=indx(node,3)) Then Write (1, ) node_x(node), node_y(node) End If End Do Close (Unit=1) Write (, '(a)') ' ' Write (, '(a)') 'XY3_WRITE:' Write (, '(a)') ' The 3-node triangle nodal coordinates were' Write (, '(a)') ' written to the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') file_name Return 10 Continue Write (, '(a)') ' ' Write (, '(a)') 'XY3_WRITE - Fatal errorc' Write (, '(a)') ' Could not open the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') file_name Stop End Subroutine xy3_write Subroutine xy6_write(file_name, node_num, node_x, node_y) !********************************************************************72 ! !c XY6_WRITE writes the coordinates of nodes of 6-node triangles to a file. ! Implicit None Integer node_num Character (Len=) file_name Integer node Double Precision node_x(node_num) Double Precision node_y(node_num) Open (Unit=1, File=file_name, Err=10) Do node = 1, node_num Write (1, ) node_x(node), node_y(node) End Do Close (Unit=1) Write (, '(a)') ' ' Write (, '(a)') 'XY6_WRITE:' Write (, '(a)') ' The 6-node triangle nodal coordinates were' Write (, '(a)') ' written to the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') file_name Return 10 Continue Write (, '(a)') ' ' Write (, '(a)') 'XY6_WRITE - Fatal errorc' Write (, '(a)') ' Could not open the output file:' Write (, '(a)') ' ' Write (, '(4x,a)') file_name Stop End Subroutine xy6_write !c$$$Run, Output: HCELL: FORTRAN95,2013&18 version Solve the Navier Stokes fluid flow equations in an H-shaped region, using finite elements. The X direction is divided into three regions, with element densities: `45 15 45` Corresponding NX = 106 The Y direction is divided into three regions, with element densities: `5 1 5` Corresponding NY = 12 The X subregions are demarked by 4 values: `0.0000 45.0000 60.0000 105.0000` The Y subregions are demarked by 4 values: `0.0000 5.0000 6.0000 11.0000` Maximum number of nodes = 4853 Maximum number of elements = 2310 Maximum number of unknowns = 10978 HCELL_NODE_COUNT: Number of nodes = 4673 HCELL_ELEMENT_COUNT: Number of elements = 2130 HCELL_DOF_COUNT: Number of degrees of freedom = 10618 ELEMENT_NODE_BANDWIDTH: Lower half bandwidth = 106 Upper half bandwidth = 106 Total bandwidth = 213 Number of nodes = 4673 Number of elements = 2130 Number of unknowns = 10618 XY6_WRITE: The 6-node triangle nodal coordinates were written to the output file: `xy6.txt` XY3_WRITE: The 3-node triangle nodal coordinates were written to the output file: `xy3.txt` ELEMENT_NODE_WRITE: The 3-node triangle nodal coordinates were written to the output file: `element_node.txt` NODE_EPS: An encapsulated Post file was created containing an image of the nodes. The file name is hcell_nodes.eps ...Program finished with exit code 0 Press ENTER to exit console. 0
2022-11-25 10:46:03　　评论淘帖0 举报相关推荐 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之a 753 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b 729 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之c 294 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之图形表达 580 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之b 362 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之a 392 • 纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域 475 • 物理学最难的方程之一纳维-斯托克斯方程，解答奖金高达100万美金 27944 • 有限元法的原理 1194 • 精准估算压缩性稳态斯托克斯方程于2D单位正方形 226