CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
hyperelastics.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine hyperelastics (inpc, textpart, elcon, nelcon, nmat, ntmat_, ncmat_, irstrt, istep, istat, n, iperturb, iline, ipol, inl, ipoinp, inp, ipoinpc)
 

Function/Subroutine Documentation

◆ hyperelastics()

subroutine hyperelastics ( character*1, dimension(*)  inpc,
character*132, dimension(16)  textpart,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
integer  nmat,
integer  ntmat_,
integer  ncmat_,
integer  irstrt,
integer  istep,
integer  istat,
integer  n,
integer, dimension(*)  iperturb,
integer  iline,
integer  ipol,
integer  inl,
integer, dimension(2,*)  ipoinp,
integer, dimension(3,*)  inp,
integer, dimension(0:*)  ipoinpc 
)
22 !
23 ! reading the input deck: *HYPERELASTIC
24 !
25  implicit none
26 !
27  character*1 inpc(*)
28  character*132 textpart(16)
29 !
30  integer nelcon(2,*),nmat,ntmat,ntmat_,istep,istat,ipoinpc(0:*),
31  & n,key,i,j,k,ityp,iperturb(*),iend,jcoef(3,14),ncmat_,irstrt,
32  & iline,ipol,inl,ipoinp(2,*),inp(3,*)
33 !
34  real*8 elcon(0:ncmat_,ntmat_,*),um
35 !
36 ! jcoef indicates for each hyperelastic model the position of
37 ! the compressibility coefficients in the field elcon (max. 3
38 ! positions per model)
39 !
40  data jcoef /3,0,0,3,0,0,2,0,0,3,0,0,5,6,0,7,8,9,3,0,0,
41  & 6,7,0,12,13,14,2,0,0,3,4,0,4,5,6,5,0,0,4,5,6/
42 !
43  ntmat=0
44  iperturb(1)=3
45  iperturb(2)=1
46  write(*,*) '*INFO reading *HYPERELASTIC: nonlinear geometric'
47  write(*,*) ' effects are turned on'
48  write(*,*)
49 !
50  if((istep.gt.0).and.(irstrt.ge.0)) then
51  write(*,*)
52  & '*ERROR reading *HYPERELASTIC: *HYPERELASTIC should be'
53  write(*,*) ' placed before all step definitions'
54  call exit(201)
55  endif
56 !
57  if(nmat.eq.0) then
58  write(*,*)
59  & '*ERROR reading *HYPERELASTIC: *HYPERELASTIC should be'
60  write(*,*) ' preceded by a *MATERIAL card'
61  call exit(201)
62  endif
63 !
64  ityp=-7
65 !
66  do i=2,n
67  if(textpart(i)(1:12).eq.'ARRUDA-BOYCE') then
68  ityp=-1
69  elseif(textpart(i)(1:13).eq.'MOONEY-RIVLIN') then
70  ityp=-2
71  elseif(textpart(i)(1:8).eq.'NEOHOOKE') then
72  ityp=-3
73  elseif(textpart(i)(1:5).eq.'OGDEN') then
74  ityp=-4
75  elseif(textpart(i)(1:10).eq.'POLYNOMIAL') then
76  ityp=-7
77  elseif(textpart(i)(1:17).eq.'REDUCEDPOLYNOMIAL') then
78  ityp=-10
79  elseif(textpart(i)(1:11).eq.'VANDERWAALS') then
80  ityp=-13
81  elseif(textpart(i)(1:4).eq.'YEOH') then
82  ityp=-14
83  elseif(textpart(i)(1:2).eq.'N=') then
84  if(textpart(i)(3:3).eq.'1') then
85  elseif(textpart(i)(3:3).eq.'2') then
86  if(ityp.eq.-4) then
87  ityp=-5
88  elseif(ityp.eq.-7) then
89  ityp=-8
90  elseif(ityp.eq.-10) then
91  ityp=-11
92  else
93  write(*,*) '*WARNING reading *HYPERELASTIC: N=2 is not
94  & applicable for this material type; '
95  call inputerror(inpc,ipoinpc,iline,
96  &"*HYPERELASTIC%")
97  endif
98  elseif(textpart(i)(3:3).eq.'3') then
99  if(ityp.eq.-4) then
100  ityp=-6
101  elseif(ityp.eq.-7) then
102  ityp=-9
103  elseif(ityp.eq.-10) then
104  ityp=-12
105  else
106  write(*,*) '*WARNING reading *HYPERELASTIC: N=3 is not
107  & applicable for this material type; '
108  call inputerror(inpc,ipoinpc,iline,
109  &"*HYPERELASTIC%")
110  endif
111  else
112  write(*,*) '*WARNING reading *HYPERELASTIC: only N=1, N=2
113  &, or N=3 are allowed; '
114  call inputerror(inpc,ipoinpc,iline,
115  &"*HYPERELASTIC%")
116  endif
117  else
118  write(*,*)
119  & '*WARNING reading *HYPERELASTIC: parameter not recognized:'
120  write(*,*) ' ',
121  & textpart(i)(1:index(textpart(i),' ')-1)
122  call inputwarning(inpc,ipoinpc,iline,
123  &"*HYPERELASTIC%")
124  endif
125  enddo
126 !
127  nelcon(1,nmat)=ityp
128 !
129  if((ityp.ne.-6).and.(ityp.ne.-9)) then
130  if((ityp.eq.-3).or.(ityp.eq.-10)) then
131  iend=2
132  elseif((ityp.eq.-1).or.(ityp.eq.-2).or.(ityp.eq.-4).or.
133  & (ityp.eq.-7)) then
134  iend=3
135  elseif(ityp.eq.-11) then
136  iend=4
137  elseif(ityp.eq.-13) then
138  iend=5
139  elseif((ityp.eq.-5).or.(ityp.eq.-12).or.(ityp.eq.-14)) then
140  iend=6
141  elseif(ityp.eq.-8) then
142  iend=7
143  endif
144  do
145  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
146  & ipoinp,inp,ipoinpc)
147  if((istat.lt.0).or.(key.eq.1)) exit
148  ntmat=ntmat+1
149  nelcon(2,nmat)=ntmat
150  if(ntmat.gt.ntmat_) then
151  write(*,*)'*ERROR reading *HYPERELASTIC: increase ntmat_'
152  call exit(201)
153  endif
154  do i=1,iend
155  read(textpart(i)(1:20),'(f20.0)',iostat=istat)
156  & elcon(i,ntmat,nmat)
157  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
158  &"*HYPERELASTIC%")
159  enddo
160  read(textpart(iend+1)(1:20),'(f20.0)',iostat=istat)
161  & elcon(0,ntmat,nmat)
162  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
163  &"*HYPERELASTIC%")
164  enddo
165  else
166  do
167  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
168  & ipoinp,inp,ipoinpc)
169  if((istat.lt.0).or.(key.eq.1)) exit
170  ntmat=ntmat+1
171  nelcon(2,nmat)=ntmat
172  if(ntmat.gt.ntmat_) then
173  write(*,*)'*ERROR reading *HYPERELASTIC: increase ntmat_'
174  call exit(201)
175  endif
176  do i=1,8
177  read(textpart(i)(1:20),'(f20.0)',iostat=istat)
178  & elcon(i,ntmat,nmat)
179  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
180  &"*HYPERELASTIC%")
181  enddo
182 !
183  if(ityp.eq.-6) then
184  iend=1
185  else
186  iend=4
187  endif
188  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
189  & ipoinp,inp,ipoinpc)
190  if((istat.lt.0).or.(key.eq.1)) then
191  write(*,*)
192  & '*ERROR reading *HYPERELASTIC: hyperelastic definition'
193  write(*,*) ' is not complete. '
194  call inputerror(inpc,ipoinpc,iline,
195  &"*HYPERELASTIC%")
196  call exit(201)
197  endif
198  do i=1,iend
199  read(textpart(i)(1:20),'(f20.0)',iostat=istat)
200  & elcon(8+i,ntmat,nmat)
201  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
202  &"*HYPERELASTIC%")
203  enddo
204  read(textpart(iend+1)(1:20),'(f20.0)',iostat=istat)
205  & elcon(0,ntmat,nmat)
206  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
207  &"*HYPERELASTIC%")
208  enddo
209  endif
210 !
211 ! if any of the compressibility coefficients is zero (incompressible
212 ! material), it is replaced. The lowest order coefficient is replaced
213 ! such that it corresponds to a Poisson coeffient of 0.475, the
214 ! following ones are replaced by a power of the first one
215 !
216  do j=1,ntmat
217 !
218 ! calculating the shear coefficient in the undeformed state
219 !
220  if(ityp.eq.-1) then
221  um=elcon(1,j,nmat)
222  elseif(ityp.eq.-2) then
223  um=2.d0*(elcon(1,j,nmat)+elcon(2,j,nmat))
224  elseif(ityp.eq.-3) then
225  um=2.d0*elcon(1,j,nmat)
226  elseif(ityp.eq.-4) then
227  um=elcon(1,j,nmat)
228  elseif(ityp.eq.-5) then
229  um=elcon(1,j,nmat)+elcon(3,j,nmat)
230  elseif(ityp.eq.-6) then
231  um=elcon(1,j,nmat)+elcon(3,j,nmat)+elcon(5,j,nmat)
232  elseif((ityp.eq.-7).or.(ityp.eq.-8).or.(ityp.eq.-9)) then
233  um=2.d0*(elcon(1,j,nmat)+elcon(2,j,nmat))
234  elseif((ityp.eq.-10).or.(ityp.eq.-11).or.(ityp.eq.-12)) then
235  um=2.d0*elcon(1,j,nmat)
236  elseif(ityp.eq.-13) then
237  um=elcon(1,j,nmat)
238  elseif(ityp.eq.-14) then
239  um=2.d0*elcon(1,j,nmat)
240  endif
241 !
242  do i=1,3
243  k=jcoef(i,abs(ityp))
244  if(k.eq.0) exit
245  if(dabs(elcon(k,j,nmat)).lt.1.d-10) then
246  elcon(k,j,nmat)=(0.1d0/um)**i
247  write(*,*)
248  & '*WARNING reading *HYPERELASTIC: default value was'
249  write(*,*) ' used for compressibility coefficient
250  &s'
251  write(*,100) i,elcon(k,j,nmat)
252  endif
253  enddo
254  enddo
255 !
256  100 format(' D',i1,' = ',e11.4)
257 !
258  return
subroutine inputwarning(inpc, ipoinpc, iline, text)
Definition: inputwarning.f:20
subroutine getnewline(inpc, textpart, istat, n, key, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: getnewline.f:21
static double * e11
Definition: radflowload.c:42
subroutine inputerror(inpc, ipoinpc, iline, text)
Definition: inputerror.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)