MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
foamgeom.f90
Go to the documentation of this file.
1 
7 module foamgeom
8  use constants
9  use solve_nonlin
10  implicit none
11  private
12  public foam_morpholgy
13 contains
14 !********************************BEGINNING*************************************
16 subroutine foam_morpholgy
17  integer, parameter :: n=2
18  integer :: info,i
19  real (dp) :: tol=1e-8_dp
20  real (dp), dimension(n) :: x,fvec,diag
21  write(*,*) 'Foam morphology:'
22  write(mfi,*) 'Foam morphology:'
23  select case(morph_input)
24  case(1) !dwall is input; fs and dstrut are calculated
25  x(1)=fs
26  x(2)=dstrut
27  call hbrd(fcn_dwall,n,x,fvec,epsilon(pi),tol,info,diag)
28  if (info /= 1) then
29  write(*,*) 'unable to determine foam morphology parameters, &
30  hbrd returned',info
31  write(mfi,*) 'unable to determine foam morphology parameters, &
32  hbrd returned',info
33  stop
34  endif
35  fs=x(1)
36  dstrut=x(2)
37  if (fs<0 .or. dstrut< 0) then
38  write(*,*) 'unable to determine foam morphology &
39  parameters, try different initial guess'
40  write(mfi,*) 'unable to determine foam morphology &
41  parameters, try different initial guess'
42  stop
43  endif
44  case(2) !fs is input; dwall and dstrut are calculated
45  if (fs<struttol) then
46  dwall=(1-por)*dcell/3.775_dp
47  dstrut=0
48  else
49  do i=1,10
50  x(1)=dwall*i
51  x(2)=dstrut*i
52  call hbrd(fcn_fs,n,x,fvec,epsilon(pi),tol,info,diag)
53  if (info /= 1) then
54  write(*,*) 'unable to determine foam morphology &
55  parameters, hbrd returned',info, 'restarting'
56  write(mfi,*) 'unable to determine foam morphology &
57  parameters, hbrd returned',info, 'restarting'
58  else
59  exit
60  endif
61  enddo
62  if (info /= 1) then
63  write(*,*) 'unable to determine foam morphology parameters, &
64  hbrd returned',info
65  write(mfi,*) 'unable to determine foam morphology parameters, &
66  hbrd returned',info
67  stop
68  endif
69  dwall=x(1)
70  dstrut=x(2)
71  if (dwall<0 .or. dstrut< 0) then
72  write(*,*) 'unable to determine foam &
73  morphology parameters, try different initial guess'
74  write(mfi,*) 'unable to determine foam &
75  morphology parameters, try different initial guess'
76  stop
77  endif
78  endif
79  case(3) !dstrut is input; dwall and fs are calculated
80  x(1)=dwall
81  x(2)=fs
82  call hbrd(fcn_dstrut,n,x,fvec,epsilon(pi),tol,info,diag)
83  if (info /= 1) then
84  write(*,*) 'unable to determine foam morphology parameters, &
85  hbrd returned',info
86  write(mfi,*) 'unable to determine foam morphology parameters, &
87  hbrd returned',info
88  stop
89  endif
90  dwall=x(1)
91  fs=x(2)
92  if (dwall<0 .or. fs< 0) then
93  write(*,*) 'unable to determine foam morphology &
94  parameters, try different initial guess'
95  write(mfi,*) 'unable to determine foam morphology &
96  parameters, try different initial guess'
97  stop
98  endif
99  case(4) !fs is input; dwall and dstrut are calculated
100  if (fs<struttol) then
101  dwall=(1-por)*dcell/3.775_dp
102  dstrut=0
103  else
104  x(1)=dwall
105  x(2)=dstrut
106  call hbrd(fcn_fs2,n,x,fvec,epsilon(pi),tol,info,diag)
107  if (info /= 1) then
108  write(*,*) 'unable to determine foam morphology parameters, &
109  hbrd returned',info
110  write(mfi,*) 'unable to determine foam morphology parameters, &
111  hbrd returned',info
112  stop
113  endif
114  dwall=x(1)
115  dstrut=x(2)
116  if (dwall<0 .or. dstrut< 0) then
117  write(*,*) 'unable to determine foam &
118  morphology parameters, try different initial guess'
119  write(mfi,*) 'unable to determine foam &
120  morphology parameters, try different initial guess'
121  stop
122  endif
123  endif
124  case default
125  write(*,*) 'unknown foam morphology input'
126  write(mfi,*) 'unknown foam morphology input'
127  stop
128  end select
129  write(*,'(2x,A,1x,e9.3)') 'porosity:', por
130  write(*,'(2x,A,1x,e9.3,1x,A)') 'foam density:', rhof, 'kg/m^3'
131  write(*,'(2x,A,1x,e9.3,1x,A)') 'cell size:', dcell*1e6, 'um'
132  write(*,'(2x,A,1x,e9.3,1x,A)') 'wall thickness:', dwall*1e6, 'um'
133  write(*,'(2x,A,1x,e9.3)') 'strut content:', fs
134  write(*,'(2x,A,1x,e9.3,1x,A)') 'strut diameter:', dstrut*1e6, 'um'
135  write(mfi,'(2x,A,1x,e9.3)') 'porosity:', por
136  write(mfi,'(2x,A,1x,e9.3,1x,A)') 'foam density:', rhof, 'kg/m^3'
137  write(mfi,'(2x,A,1x,e9.3,1x,A)') 'cell size:', dcell*1e6, 'um'
138  write(mfi,'(2x,A,1x,e9.3,1x,A)') 'wall thickness:', dwall*1e6, 'um'
139  write(mfi,'(2x,A,1x,e9.3)') 'strut content:', fs
140  write(mfi,'(2x,A,1x,e9.3,1x,A)') 'strut diameter:', dstrut*1e6, 'um'
141 end subroutine foam_morpholgy
142 !***********************************END****************************************
143 
144 
145 !********************************BEGINNING*************************************
147 subroutine fcn_dwall(n,x,fvec,iflag)
148  integer, intent(in) :: n
149  real (dp), intent(in) :: x(n)
150  real (dp), intent(out) :: fvec(n)
151  integer, intent(inout) :: iflag
152  real(dp) :: Vcell,Vstruts,Vwalls,fs,dstrut,dcelldd
153  fs=x(1)
154  dstrut=x(2)
155  dcelldd=dcell*(pi/6/0.348_dp)**(1/3._dp)
156  vcell=0.348_dp*dcelldd**3
157  vstruts=2.8_dp*dstrut**2*dcelldd-3.93_dp*dstrut**3
158  vwalls=(1.3143_dp*dcelldd**2-7.367_dp*dstrut*dcelldd+10.323_dp*dstrut**2)*&
159  dwall
160  fvec(1)=fs-vstruts/(vstruts+vwalls)
161  fvec(2)=1-por-(vstruts+vwalls)/vcell
162 end subroutine fcn_dwall
163 !***********************************END****************************************
164 
165 
166 !********************************BEGINNING*************************************
168 subroutine fcn_fs(n,x,fvec,iflag)
169  integer, intent(in) :: n
170  real (dp), intent(in) :: x(n)
171  real (dp), intent(out) :: fvec(n)
172  integer, intent(inout) :: iflag
173  real(dp) :: Vcell,Vstruts,Vwalls,dwall,dstrut,dcelldd
174  dwall=x(1)
175  dstrut=x(2)
176  dcelldd=dcell*(pi/6/0.348_dp)**(1/3._dp)
177  vcell=0.348_dp*dcelldd**3
178  vstruts=2.8_dp*dstrut**2*dcelldd-3.93_dp*dstrut**3
179  vwalls=(1.3143_dp*dcelldd**2-7.367_dp*dstrut*dcelldd+10.323_dp*dstrut**2)*&
180  dwall
181  fvec(1)=fs-vstruts/(vstruts+vwalls)
182  fvec(2)=1-por-(vstruts+vwalls)/vcell
183 end subroutine fcn_fs
184 !***********************************END****************************************
185 
186 
187 !********************************BEGINNING*************************************
189 subroutine fcn_dstrut(n,x,fvec,iflag)
190  integer, intent(in) :: n
191  real (dp), intent(in) :: x(n)
192  real (dp), intent(out) :: fvec(n)
193  integer, intent(inout) :: iflag
194  real(dp) :: Vcell,Vstruts,Vwalls,dwall,fs,dcelldd
195  dwall=x(1)
196  fs=x(2)
197  dcelldd=dcell*(pi/6/0.348_dp)**(1/3._dp)
198  vcell=0.348_dp*dcelldd**3
199  vstruts=2.8_dp*dstrut**2*dcelldd-3.93_dp*dstrut**3
200  vwalls=(1.3143_dp*dcelldd**2-7.367_dp*dstrut*dcelldd+10.323_dp*dstrut**2)*&
201  dwall
202  fvec(1)=fs-vstruts/(vstruts+vwalls)
203  fvec(2)=1-por-(vstruts+vwalls)/vcell
204 end subroutine fcn_dstrut
205 !***********************************END****************************************
206 
207 
208 !********************************BEGINNING*************************************
212 subroutine fcn_fs2(n,x,fvec,iflag)
213  integer, intent(in) :: n
214  real (dp), intent(in) :: x(n)
215  real (dp), intent(out) :: fvec(n)
216  integer, intent(inout) :: iflag
217  real(dp) :: Vcell,Vstruts,Vwalls,dwall,dstrut,dcelldd,x1,x2,x3
218  dwall=x(1)
219  dstrut=x(2)
220  dcelldd=dcell*(pi/6/0.348_dp)**(1/3._dp)
221  vcell=0.348_dp*dcelldd**3
222  vstruts=2.8_dp*dstrut**2*dcelldd
223  vwalls=(1.317_dp*dcelldd**2-13.4284_dp*dstrut*dcelldd+34.2375_dp*dstrut**2)*&
224  dwall+(4.639_dp*dcell-17.976_dp*dstrut)*dwall**2
225  fvec(1)=fs-vstruts/(vstruts+vwalls)
226  fvec(2)=1-por-(vstruts+vwalls)/vcell
227 end subroutine fcn_fs2
228 !***********************************END****************************************
229 end module foamgeom