teratail header banner
teratail header banner
質問するログイン新規登録

質問編集履歴

3

画像を追加

2020/04/23 13:25

投稿

highkazu
highkazu

スコア11

title CHANGED
File without changes
body CHANGED
@@ -185,4 +185,8 @@
185
185
 
186
186
  2020/04/22 追記
187
187
  法線を解析的に求めた場合の画像.
188
- ![イメージ説明](021e83ca8d57de0e6681cf8078c3a8d0.png)
188
+ ![イメージ説明](021e83ca8d57de0e6681cf8078c3a8d0.png)
189
+
190
+ 2020/04/23 追記
191
+ 正規化方法をコメント頂いたように直した場合
192
+ ![イメージ説明](ffe8aaff7f8045be5be4a085a34ceae6.png)

2

プログラムの更新

2020/04/23 13:25

投稿

highkazu
highkazu

スコア11

title CHANGED
File without changes
body CHANGED
@@ -14,33 +14,38 @@
14
14
 
15
15
  ```fortran
16
16
 
17
+
17
18
  program main
18
19
  ! Determine some arrays
19
20
  implicit none
20
- integer,parameter :: i_max=50, j_max=50,k_max=50
21
+ integer,parameter :: i_max=100, j_max=100
21
22
  double precision :: x(0:i_max,0:j_max), y(0:i_max,0:j_max), z(0:i_max,0:j_max) ! this is the position
22
23
  double precision :: v(3) ! this is the velocity (uniform flow)
23
- double precision :: nx(0:i_max-1,0:j_max-1,0:k_max-1),ny(0:i_max-1,0:j_max-1,0:k_max-1),nz(0:i_max-1,0:j_max-1,0:k_max-1)
24
24
  double precision :: cp(0:i_max,0:j_max)
25
25
  integer :: i,j
26
26
 
27
27
  call grid
28
28
  call normal_and_Cp
29
29
  call dataSave
30
+ print *, pi_generator()
30
31
  stop
31
32
  contains
32
33
 
34
+ double precision function pi_generator()
35
+ pi_generator = 4*atan(1.0)
36
+ return
37
+ end function
38
+
33
39
  subroutine grid
34
40
  ! Generating grids
35
41
  double precision :: r = 2 ! radius
36
42
  double precision :: theta(0:i_max)
37
43
  double precision :: phi(0:j_max)
38
44
  do 1000 j=0,j_max
39
- theta(j)=float(j)/float((j_max))*(3.141592)/2.0
45
+ theta(j)=float(j)/float((j_max))*(pi_generator() )/2.0
40
- print *, theta(j), j
41
46
  1000 continue
42
47
  do 1001 i=0,i_max
43
- phi(i) = float(i)/float(i_max)*3.141592*2
48
+ phi(i) = float(i)/float(i_max)*(pi_generator() )*2
44
49
  1001 continue
45
50
  ! 球座標系で,r方向の変化はないモデルなので,2変数で書けそうだ
46
51
  do 100 i=0,i_max
@@ -49,6 +54,8 @@
49
54
  y(i,j) = r*sin(phi(i))*sin(theta(j))
50
55
  z(i,j) = r*cos(theta(j))
51
56
  101 continue
57
+ ! x(i,0)=x(i,j_max)
58
+ ! y(i,0)=y(i,j_max)
52
59
  100 continue
53
60
  end subroutine
54
61
 
@@ -60,13 +67,13 @@
60
67
  ! 面と格子がずれることになるが,そこは一度妥協する
61
68
 
62
69
  !
63
- ! o (U)
70
+ ! o(U)
64
71
  ! |
72
+ ! (L) |
73
+ ! o-------o-------o(R)
65
74
  ! |
66
- !(L)o-------o-------o(R)
67
75
  ! |
68
- ! |
69
- ! o (B)
76
+ ! o(B)
70
77
 
71
78
  double precision :: vec_R(3),vec_L(3), vec_U(3), vec_B(3), &
72
79
  & cross_RU(3), cross_UL(3), cross_LB(3), cross_BR(3), &
@@ -94,7 +101,6 @@
94
101
  vec_B(1) = x(i,j-1) - x(i,j)
95
102
  vec_B(2) = y(i,j-1) - y(i,j)
96
103
  vec_B(3) = z(i,j-1) - z(i,j)
97
- print *, vec_B(2)
98
104
 
99
105
 
100
106
  ! print *, vec_B(3)
@@ -125,10 +131,16 @@
125
131
  cross_avg(2) = cross_avg(2)/sqrt(cross_avg(2)*cross_avg(2)+1.0d-5)
126
132
  cross_avg(3) = cross_avg(3)/sqrt(cross_avg(3)*cross_avg(3)+1.0d-5)
127
133
 
134
+ ! 球の場合,normal=(x,y,z)でした.x**2+y**2+z**2=r**2より,f=x**2+y**2+z**2-r**2としてgradとればokでした...
135
+ ! =====================================================================================
136
+ cross_avg(1)=0.5*x(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
137
+ cross_avg(2)=0.5*y(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
138
+ cross_avg(3)=0.5*z(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
139
+ ! ====================================================================================
128
140
  v(1)=0
129
141
  v(2)=0
130
142
  v(3)=-100*7
131
- if (dot_product(v,cross_avg) > 0) then
143
+ if (dot_product(v,cross_avg) < 0) then
132
144
  cp(i,j) = 2*((dot_product(v,cross_avg)) &
133
145
  & /sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)))**2
134
146
  else
@@ -147,7 +159,7 @@
147
159
  subroutine dataSave
148
160
  character (40):: filename
149
161
  integer :: fo2=12
150
- write(filename, "(a,i5.5,a)") "Newtonian" ! Dataファイル
162
+ write(filename, "(a,i5.5,a)") "vtkdata/Newtonian"
151
163
  open(fo2,file=filename)
152
164
  do 700 i=1,i_max
153
165
  do 701 j=1,j_max

1

頂いた回答をもとに計算しなおした画像を掲載

2020/04/22 10:03

投稿

highkazu
highkazu

スコア11

title CHANGED
File without changes
body CHANGED
@@ -169,4 +169,8 @@
169
169
  plotはgnuplotを使用して
170
170
  splot 'Newtonian' u 1:2:3:4 w pm3d
171
171
  で書いています.
172
- fortran90を使用しています
172
+ fortran90を使用しています
173
+
174
+ 2020/04/22 追記
175
+ 法線を解析的に求めた場合の画像.
176
+ ![イメージ説明](021e83ca8d57de0e6681cf8078c3a8d0.png)