[讨论]

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

2022-11-25 10:43:53 774 控制程序流体力学

0 !*********************************************************** !********************************************************* ! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b !********************************************************* ! SINOMACH ! ! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之b ! ! 周勇 20221125 !********************************************************** !c $$$ Run, Output: !********************************************************** !********************************************************** Subroutine hcell_dof_set(region_density_x, region_density_y, maxnd,& node_num, node_dof_index) !******************************************************************72 ! !c HCELL_DOF_SET assigns degrees of freedom to each node. ! Implicit None Integer maxnd Integer col Integer dof Logical found Integer i Integer j Integer node Integer node_dof_index(maxnd, 3) Integer node_num Integer region_density_x(3) Integer region_density_y(3) Integer row node = 0 dof = 0 j = 0 ! ! Working in column 1, and NOT handling nodes in extreme right. ! Do col = 1, 2region_density_x(1) ! ! Working in row 1. ! j = j + 1 i = 0 `Do row = 1, 2region_density_y(1) + 1 i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` ! ! Working in row 3. ! i = i + 2region_density_y(2) - 1 `Do row = 1, 2region_density_y(3) + 1 i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` End Do ! ! Working in column 2, including nodes on extreme right and left. ! Do col = 1, 2region_density_x(2) + 1 ! ! Working in row 1. ! j = j + 1 i = 0 `Do row = 1, 2region_density_y(1) + 1 i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` ! ! Working in row 2. ! Do row = 2, 2region_density_y(2) `i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` ! ! Working in row 3. ! Do row = 1, 2region_density_y(3) + 1 `i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` End Do ! ! Working in column 3, and NOT handling nodes in extreme left. ! Do col = 2, 2region_density_x(3) + 1 ! ! Working in row 1. ! j = j + 1 i = 0 `Do row = 1, 2region_density_y(1) + 1 i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` ! ! Working in row 3. ! i = i + 2region_density_y(2) - 1 `Do row = 1, 2region_density_y(3) + 1 i = i + 1 node = node + 1 dof = dof + 1 node_dof_index(node, 1) = dof dof = dof + 1 node_dof_index(node, 2) = dof If (mod(j,2)==1 .And. mod(i,2)==1) Then dof = dof + 1 node_dof_index(node, 3) = dof Else node_dof_index(node, 3) = 0 End If End Do` End Do ! ! Replace one continuity equation by a pressure = 0 equation. ! found = .False. Do node = 1, node_num `If (0<node_dof_index(node,3)) Then node_dof_index(node, 3) = -1 found = .True. Exit End If` End Do If (.Not. found) Then Write (, '(a)') ' ' Write (, '(a)') 'HCELL_DOF_SET - Fatal errorc' Write (, '(a)') ' Could not find a degree of freedom' Write (, '(a)') ' associated with the continuity equation.' Stop End If Return End Subroutine hcell_dof_set Subroutine hcell_element_count(region_density_x, region_density_y, element_num) !*******************************************************************72 ! !c HCELL_ELEMENT_COUNT determines the number of elements in the region. ! Implicit None Integer element_num Integer region_density_x(3) Integer region_density_y(3) element_num = 2region_density_x(1)region_density_y(1) + & 2region_density_x(1)region_density_y(3) + & 2region_density_x(2)region_density_y(1) +& 2region_density_x(2)region_density_y(2) + & 2region_density_x(2)region_density_y(3) +& 2region_density_x(3)region_density_y(1) + & 2region_density_x(3)region_density_y(3) Write (, '(a)') ' ' Write (, '(a)') 'HCELL_ELEMENT_COUNT:' Write (, '(a,i6)') ' Number of elements = ', element_num Return End Subroutine hcell_element_count Subroutine hcell_element_node(region_density_x, region_density_y,& maxel, element_num, nnodes, element_node) !********************************************************************72 ! !c HCELL_ELEMENT_NODE determines the nodes that make up each element. ! Implicit None Integer maxel Integer nnodes Integer col Integer col2 Integer element Integer element_node(maxel, nnodes) Integer element_num Integer inc1 Integer inc2 Integer node_sw Integer region_density_x(3) Integer region_density_y(3) Integer row Integer row2 element = 0 Do col = 1, 3 Do col2 = 1, region_density_x(col) Do row = 1, 3 If (row/=2 .Or. col==2) Then If (col==1) Then If (col2<region_density_x(1)) Then If (row==1) Then If (col2==1) Then node_sw = 1 Else node_sw = node_sw + inc1 + 1 End If Else node_sw = node_sw + 1 End If inc1 = (2region_density_y(1)+1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(3)+1) Else If (row==1) Then node_sw = node_sw + inc1 + 1 inc1 = (2region_density_y(1)+1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(3)+1) Else If (row==3) Then node_sw = node_sw + 1 inc1 = (2region_density_y(1)+1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(2)-1)& + (2region_density_y(3)+1) End If Else If (col==2) Then If (row==1) Then node_sw = node_sw + inc1 + 1 End If inc1 = (2region_density_y(1)+1) + (2region_density_y(2)-1)& + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(2)-1)& + (2region_density_y(3)+1) Else If (col==3) Then If (1<col2) Then If (row==1) Then node_sw = node_sw + inc1 + 1 Else node_sw = node_sw + 1 End If inc1 = (2region_density_y(1)+1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(3)+1) Else If (row==1) Then node_sw = node_sw + inc1 + 1 inc1 = (2region_density_y(1)+1) +& (2region_density_y(2)-1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(3)+1) Else If (row==3) Then node_sw = node_sw + (2region_density_y(2)-1) + 1 inc1 = (2region_density_y(1)+1) + (2region_density_y(3)+1) inc2 = (2region_density_y(1)+1) + (2region_density_y(3)+1) End If End If Do row2 = 1, region_density_y(row) element = element + 1 element_node(element, 1) = node_sw element_node(element, 2) = node_sw + inc1 + inc2 element_node(element, 3) = node_sw + 2 element_node(element, 4) = node_sw + inc1 element_node(element, 5) = node_sw + inc1 + 1 element_node(element, 6) = node_sw + 1 element = element + 1 element_node(element, 1) = node_sw + inc1 + inc2 + 2 element_node(element, 2) = node_sw + 2 element_node(element, 3) = node_sw + inc1 + inc2 element_node(element, 4) = node_sw + inc1 + 2 element_node(element, 5) = node_sw + inc1 + 1 element_node(element, 6) = node_sw + inc1 + inc2 + 1 node_sw = node_sw + 2 End Do End If End Do End Do End Do Return End Subroutine hcell_element_node Subroutine hcell_node_count(region_density_x, region_density_y, node_num) !********************************************************************72 ! !c HCELL_NODE_NUM determines the number of nodes in the region. ! Implicit None Integer node_num Integer region_density_x(3) Integer region_density_y(3) ! ! Count the nodes. ! node_num = (2region_density_x(1)+1)(2region_density_y(1)+1)& (2region_density_x(1)+1)(2region_density_y(3)+1) & (2region_density_x(2)-1)(2region_density_y(1)+1)& (2region_density_x(2)+1)(2region_density_y(2)-1)& (2region_density_x(2)-1)(2region_density_y(3)+1)& (2region_density_x(3)+1)(2region_density_y(1)+1) & (2region_density_x(3)+1)(2region_density_y(3)+1) Write (, '(a)') ' ' Write (, '(a)') 'HCELL_NODE_COUNT:' Write (, '(a,i6)') ' Number of nodes = ', node_num Return End Subroutine hcell_node_count Subroutine hcell_node_xy(region_density_x, region_density_y,& node_num, region_x, region_y, node_x, node_y) !*******************************************************************72 ! !c HCELL_NODE_XY assigns coordinates to each node. ! Implicit None Integer node_num Integer col Integer i Integer j Integer node Double Precision node_x(node_num) Double Precision node_y(node_num) Integer region_density_x(3) Integer region_density_y(3) Double Precision region_x(4) Double Precision region_y(4) Integer row node = 0 j = 0 ! ! Working in column 1, except for the extreme right. ! Do col = 1, 2region_density_x(1) `j = j + 1 i = 0` ! ! Working in row 1. ! Do row = 1, 2region_density_y(1) + 1 `i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(1)+1-col)region_x(1)+& dble(-1+col)region_x(2))/dble(2region_density_x(1)) node_y(node) = (dble(2region_density_y(1)+1-row)region_y(1)+& dble(-1+row)region_y(2))/dble(2region_density_y(1)) End Do` ! ! Working in row 3. ! i = i + 2region_density_y(2) - 1 `Do row = 1, 2region_density_y(3) + 1 i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(1)+1-col)region_x(1)+& dble(-1+col)region_x(2))/dble(2region_density_x(1)) node_y(node) = (dble(2region_density_y(3)+1-row)region_y(3)+& dble(-1+row)region_y(4))/dble(2region_density_y(3)) End Do` End Do ! ! Working in column 2, including extreme left and right. ! Do col = 1, 2region_density_x(2) + 1 ! ! Working in row 1. ! j = j + 1 i = 0 `Do row = 1, 2region_density_y(1) + 1 i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(2)+1-col)region_x(2)+& dble(-1+col)region_x(3))/dble(2region_density_x(2)) node_y(node) = (dble(2region_density_y(1)+1-row)region_y(1)+& dble(-1+row)region_y(2))/dble(2region_density_y(1)) End Do` ! ! Working in row 2. ! Do row = 2, 2region_density_y(2) `i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(2)+1-col)region_x(2)+& dble(-1+col)region_x(3))/dble(2region_density_x(2)) node_y(node) = (dble(2region_density_y(2)+1-row)region_y(2)+& dble(-1+row)region_y(3))/dble(2region_density_y(2)) End Do` ! ! Working in row 3. ! Do row = 1, 2region_density_y(3) + 1 `i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(2)+1-col)region_x(2)+& dble(-1+col)region_x(3))/dble(2region_density_x(2)) node_y(node) = (dble(2region_density_y(3)+1-row)region_y(3)+& dble(-1+row)region_y(4))/dble(2region_density_y(3)) End Do` End Do ! ! Working in column 3, except for extreme left. ! Do col = 2, 2region_density_x(3) + 1 ! ! Working in row 1. ! j = j + 1 i = 0 `Do row = 1, 2region_density_y(1) + 1 i = i + 1 node = node + 1 node_x(node) = (dble(2region_density_x(3)+1-col)region_x(3)+& dble(-1+col)region_x(4))/dble(2region_density_x(3)) node_y(node) = (dble(2region_density_y(1)+1-row)region_y(1)+& dble(-1+row)region_y(2))/dble(2region_density_y(1)) End Do i = i + 2region_density_y(2) - 1` ! ! Working in row 3. ! Do row = 1, 2region_density_y(3) + 1 i = i + 1 node = node + 1 `node_x(node) = (dble(2region_density_x(3)+1-col)& region_x(3)+dble(-1+col)region_x(4))/dble(2region_density_x(3)) node_y(node) = (dble(2region_density_y(3)+1-row)& region_y(3)+dble(-1+row)region_y(4))/dble(2region_density_y(3)) End Do` End Do Return End Subroutine hcell_node_xy Integer Function idamax(n, dx, incx) !******************************************************************72 ! !c IDAMAX finds the vector element of largest magnitude. ! Double Precision dx(1), dmax Integer i, incx, ix, n idamax = 0 If (n<1) Return idamax = 1 If (n==1) Return If (incx==1) Goto 20 ! ! CODE FOR INCREMENT NOT EQUAL TO 1 ! ix = 1 dmax = dabs(dx(1)) ix = ix + incx Do i = 2, n If (dabs(dx(ix))<=dmax) Goto 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx End Do Return ! ! CODE FOR INCREMENT EQUAL TO 1 ! 20 dmax = dabs(dx(1)) Do i = 2, n If (dabs(dx(i))<=dmax) Goto 30 idamax = i dmax = dabs(dx(i)) 30 End Do Return End Function idamax Subroutine node_eps(node_file_name, np, node_mask, xc, yc, title) !******************************************************************72 ! !c NODE_EPS creates an EPS file containing an image of the nodes. ! Implicit None Integer np Double Precision ave_x Double Precision ave_y Integer circle_size Double Precision dif Integer eps_unit Integer eps_x Integer eps_y Character (Len=) node_file_name Integer i Integer ios Integer ip1 Integer j Integer k Integer node Logical node_mask(np) Integer node_unmasked_num Double Precision node_x_max Double Precision node_x_min Double Precision node_y_max Double Precision node_y_min Double Precision r8_huge Double Precision scale Character (Len=40) string Character (Len=) title Double Precision xc(np) Double Precision yc(np) ! ! Determine the range and number of the unmasked nodes. ! node_x_min = r8_huge() node_x_max = -r8_huge() node_y_min = r8_huge() node_y_max = -r8_huge() node_unmasked_num = 0 Do node = 1, np If (node_mask(node)) Then node_unmasked_num = node_unmasked_num + 1 node_x_min = min(node_x_min, xc(node)) node_x_max = max(node_x_max, xc(node)) node_y_min = min(node_y_min, yc(node)) node_y_max = max(node_y_max, yc(node)) End If End Do If (node_unmasked_num<=200) Then circle_size = 3 Else If (node_unmasked_num<=800) Then circle_size = 2 Else circle_size = 1 End If If (node_y_max-node_y_min<node_x_max-node_x_min) Then scale = node_x_max - node_x_min dif = (node_x_max-node_x_min) - (node_y_max-node_y_min) node_y_max = node_y_max + 0.5dif node_y_min = node_y_min - 0.5dif Else scale = node_y_max - node_y_min dif = (node_y_max-node_y_min) - (node_x_max-node_x_min) node_x_max = node_x_max + 0.5dif node_x_min = node_x_min - 0.5dif End If eps_unit = 1 Open (Unit=eps_unit, File=node_file_name, Status='unknown') Write (eps_unit, '(a)') '%cPS-Adobe-3.0 EPSF-3.0' Write (eps_unit, '(a)') '%%Creator: node_eps(hcell.f)' Write (eps_unit, '(a,a)') '%% : ', node_file_name Write (eps_unit, '(a)') '%%Pages: 1' Write (eps_unit, '(a)') '%%BoundingBox: 36 36 576 756' Write (eps_unit, '(a)') '%%Document-Fonts: Times-Roman' Write (eps_unit, '(a)') '%%LanguageLevel: 1' Write (eps_unit, '(a)') '%%EndComments' Write (eps_unit, '(a)') '%%BeginProlog' Write (eps_unit, '(a)') '/inch {72 mul} def' Write (eps_unit, '(a)') '%%EndProlog' Write (eps_unit, '(a)') '%%Page: 1 1' Write (eps_unit, '(a)') 'save' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Set RGB line color.' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 0.9000 0.9000 0.9000 setrgbcolor' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Draw a gray border around the page.' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') 'newpath' Write (eps_unit, '(a)') ' 36 126 moveto' Write (eps_unit, '(a)') ' 576 126 lineto' Write (eps_unit, '(a)') ' 576 666 lineto' Write (eps_unit, '(a)') ' 36 666 lineto' Write (eps_unit, '(a)') ' 36 126 lineto' Write (eps_unit, '(a)') 'stroke' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Set RGB line color.' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Label the plot:' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor' Write (eps_unit, '(a)') '/Times-Roman findfont 0.50 inch scalefont setfont' Write (eps_unit, '(a)') ' 36 666 moveto' Write (eps_unit, '(a,a,a)') '(',title, ') show' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Define a clipping polygon' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 36 126 moveto' Write (eps_unit, '(a)') ' 576 126 lineto' Write (eps_unit, '(a)') ' 576 666 lineto' Write (eps_unit, '(a)') ' 36 666 lineto' Write (eps_unit, '(a)') ' 36 126 lineto' Write (eps_unit, '(a)') 'clip newpath' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Draw filled dots at each node:' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 0.0000 0.0000 1.0000 setrgbcolor' Do node = 1, np `If (node_mask(node)) Then eps_x = int((node_x_max-xc(node))61.0+(+xc(node)-node_x_min)551.0)/scale eps_y = int((node_y_max-yc(node))151.0+(yc(node)-node_y_min)641.0)/scale Write (eps_unit, '(a,i4,2x,i4,2x,i4,a)') 'newpath ',& eps_x, eps_y, circle_size, ' 0 360 arc closepath fill' End If` End Do ! ! Label each node, but only if there aren't too many of themc ! If (node_unmasked_num<=200) Then Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% Label the nodes:' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor' Write (eps_unit, '(a)') '/Times-Roman findfont 0.20 inch scalefont setfont' Do node = 1, np If (node_mask(node)) Then eps_x = int((node_x_max-xc(node))61.0+(+xc(node)-node_x_min)551.0)/scale eps_y = int((node_y_max-yc(node))151.0+(yc(node)-node_y_min)641.0)/scale Write (eps_unit, '(i4,2x,i4,a,i4,a)') eps_x, eps_y + 5,& ' moveto (', node, ') show' End If End Do End If Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') 'restore showpage' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '% End of page' Write (eps_unit, '(a)') '%' Write (eps_unit, '(a)') '%%Trailer' Write (eps_unit, '(a)') '%%EOF' Close (Unit=eps_unit) Write (, '(a)') ' ' Write (, '(a)') 'NODE_EPS:' Write (, '(a)') ' An encapsulated Post file was created' Write (, '(a)') ' containing an image of the nodes.' Write (, '(a)') ' The file name is' Write (, '(4x,a)') node_file_name Return End Subroutine node_eps Subroutine nstoke(xc, yc, area, xm, ym, a, f, g, uold, reynld,& tolns, xlngth, ylngth, node, indx, ipivot, mrow1, nlband, & nuband, nband, nrow1, ncol1, nelemn, np, nnodes, nuk, nquad,& neqn, nsteps, nsim, maxnd, maxel, rdel, alpha) !*******************************************************************72 ! !c NSTOKE solves the Navier-Stokes equations using Taylor-Hood elements. ! Implicit Double Precision (A-H, O-Z) Double Precision a(mrow1, ) Double Precision area() Double Precision f() Double Precision g() Integer indx(maxnd, 3) Integer ipivot() Integer it Integer iter Integer node(maxel, ) Double Precision reynld Double Precision un(2) Double Precision unx(2) Double Precision uny(2) Double Precision uold() Double Precision uold_qp Double Precision visc Double Precision vold_qp Double Precision x Double Precision xc() Double Precision xm(maxel, ) Double Precision y Double Precision yc() Double Precision ym(maxel, ) ! ! zero arrays ! g array contains the previous iterate ! f array contains the right hand side initially and then the current ! iterate is overwritten on f ! visc = 1.D0/reynld ! ! matrix assembly triangle by triangle ! nsim is the number of simple iterations performed ! csim = 0.D0 Do iter = 1, nsteps niter = iter If (iter>nsim) csim = 1.D0 Do i = 1, nrow1 Do j = 1, ncol1 a(i, j) = 0.D0 End Do End Do Do it = 1, nelemn arr = area(it)/3.D0 Do iquad = 1, nquad y = ym(it, iquad) x = xm(it, iquad) Call trans(it, x, y, det, xix, xiy, etax, etay, xc, yc, node, maxel) ar = arrdet Do kk = 1, 2 un(kk) = 0.D0 uny(kk) = 0.D0 unx(kk) = 0.D0 End Do uold_qp = 0.0D0 vold_qp = 0.0D0 Do iq = 1, nnodes Call refqbf(x, y, iq, bb, tbx, tby) bx = tbxxix + tbyetax by = tbxxiy + tbyetay ip = node(it, iq) Do iuk = 1, 2 iun = indx(ip, iuk) If (0<iun) Then un(iuk) = un(iuk) + bbg(iun) unx(iuk) = unx(iuk) + bxg(iun) uny(iuk) = uny(iuk) + byg(iun) If (iuk==1) uold_qp = uold_qp + bbuold(iun) If (iuk==2) vold_qp = vold_qp + bbuold(iun) Else If (iun<0) Then ubc = alphaubdry(iuk, ip, xc, yc) un(iuk) = un(iuk) + bbubc unx(iuk) = unx(iuk) + bxubc uny(iuk) = uny(iuk) + byubc If (iuk==1) uold_qp = uold_qp + bbubc If (iuk==2) vold_qp = vold_qp + bbubc End If End Do End Do Do iq = 1, nnodes ip = node(it, iq) Call refqbf(x, y, iq, bb, tbx, tby) bx = tbxxix + tbyetax by = tbxxiy + tbyetay bbl = refbsp(x, y, iq) Do iuk = 1, nuk i = indx(ip, iuk) If (i<=0) Goto 210 If (iuk==1) f(i) = f(i) + csim((un(1)unx(1)+un(2)uny(1))& bb)ar + rdeluold_qpbbar If (iuk==2) f(i) = f(i) + csim((un(1)unx(2)+un(2)uny(2))& bb)ar + rdelvold_qpbbar Do iqq = 1, nnodes ipp = node(it, iqq) Call refqbf(x, y, iqq, bbb, tbbx, tbby) bbx = tbbxxix + tbbyetax bby = tbbxxiy + tbbyetay bbbl = refbsp(x, y, iqq) Do iukk = 1, nuk j = indx(ipp, iukk) If (j==0) Goto 190 aij = 0.D0 If (i==neqn) Goto 190 If (iuk==1) Then If (iukk==1) aij = visc(bybby+bxbbx)& + (bbbunx(1)bb)csim + bbbbxun(1) +& bbbbyun(2) + rdel(bbbbb) If (iukk==2) aij = csim(bbbbbuny(1)) If (iukk==3) aij = -bxbbbl Else If (iuk==2) Then If (iukk==2) aij = (visc(bybby+bxbbx)+& (bbbbbuny(2))csim+bbbbyun(2)+bbbbxun(1)) + rdel(bbbbb) If (iukk==1) aij = csim(bbbbbunx(2)) If (iukk==3) aij = -bybbbl Else If (iukk==1) aij = bbxbbl If (iukk==2) aij = bbybbl End If If (j<0) Goto 185 iuse = i - j + nband a(iuse, j) = a(iuse, j) + aijar Goto 190 ! ! add terms to rhs for inhomogeneous boundary condition ! 185 Continue f(i) = f(i) - aralphaubdry(iukk, ipp, xc, yc)aij 190 End Do End Do 210 End Do End Do End Do End Do f(neqn) = 0.D0 Do j = neqn - nlband, neqn - 1 i = neqn - j + nband a(i, j) = 0.D0 End Do a(nband, neqn) = 1.D0 ! ! solve system ! job = 0 Call dgbfa(a, mrow1, neqn, nlband, nuband, ipivot, info) Call dgbsl(a, mrow1, neqn, nlband, nuband, ipivot, f, job) ! ! check for convergence ! diff = 0.D0 Do i = 1, neqn diff = diff + (g(i)-f(i))*2 End Do diff = sqrt(diff) `Write (, *) ' Iteration ', iter, ' Difference is ', diff If (diff<=tolns) Then Return End If Do i = 1, neqn g(i) = f(i) f(i) = 0.D0 End Do` End Do Return End Subroutine nstoke 0
2022-11-25 10:43:53　　评论淘帖0 举报相关推荐 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之a 793 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之c 725 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之b 405 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之图形表达 596 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之c 307 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之a 438 • 纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域 519 • 物理学最难的方程之一纳维-斯托克斯方程，解答奖金高达100万美金 28024 • 有限元法的原理 1211 • 精准估算压缩性稳态斯托克斯方程于2D单位正方形 239