1 : C CONJP2 SOURCE CHAT 05/01/12 22:17:36 5004
2 : SUBROUTINE CONJP2(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT,
3 : & ga)
4 : C************************************************************************
5 : C
6 : C PROJET : CASTEM 2000
7 : C
8 : C NOM : CONJP2
9 : C
10 : C DESCRIPTION : Voir KONJP2
11 : C
12 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
13 : C
14 : C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
15 : C
16 : C************************************************************************
17 : C
18 : c----------------------------------------------------------------------
19 : c GENERAL DESCRIPTION:
20 : c This subroutine provides the jacobians which are the derivatives
21 : c of the numerical flux function defined at the cell interface
22 : c with respect to the primitive variables of the left and right
23 : c cells (relative to the cell interface).
24 : c
25 : c EQUATIONS: 2D Euler equations of gas dynamics
26 : c
27 : c
28 : c REFERENCE: JCP, 129, 364-382 (1996)
29 : c " A Sequel to AUSM: AUSM+ ".
30 : c----------------------------------------------------------------------
31 : c INPUT:
32 : c
33 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
34 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
35 : c
36 : c beta -- parameter of the AUSM+ scheme in the Mach function;
37 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
38 : c
39 : c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
40 : c left cell;
41 : c
42 : c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
43 : c right cell;
44 : c
45 : c nvect -- normal vector to the interface (2 components in 2D);
46 : c
47 : c tvect -- tangential vector to the interface;
48 : c
49 : c ga -- ratio of the specific heats (assumed constant)
50 : c----------------------------------------------------------------------
51 : c
52 : c OUTPUT:
53 : c
54 : c jl -- jakobian matrix 4 by 4 - derivatives of the numerical
55 : c flux function with respect to the primitive variables
56 : c from the left cell;
57 : c
58 : c jr -- jakobian matrix 4 by 4 - derivatives of the numerical
59 : c flux function with respect to the primitive variables
60 : c from the right cell.
61 : c
62 : c N.B. Numerical flux function is given by
63 : c \rho un
64 : c \rho un ux + p (n.x)
65 : c \rho un uy + p (n.y)
66 : c \rho un ht
67 : c----------------------------------------------------------------------
68 : IMPLICIT INTEGER(I-N)
69 : real*8 wvec_l(4),wvec_r(4)
70 : real*8 nvect(2),tvect(2)
71 : real*8 jl(4,4),jr(4,4),f(4)
72 : real*8 alpha,beta
73 : real*8 ga, gm1,temph
74 : real*8 n_x,n_y,t_x,t_y
75 : real*8 un_l, un_r, ut_l, ut_r
76 : real*8 ml,mr,Mplus,Mmin,mmid
77 : real*8 mpl_m, mmin_m,am
78 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
79 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
80 : real*8 Pplus,Pmin,pmid
81 : real*8 hr_l,hr_r,det
82 : real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
83 : real*8 aleft, arigh
84 : real*8 damr_l,damr_r,damu_l,damu_r
85 : real*8 damv_l,damv_r,damp_l,damp_r
86 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
87 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
88 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
89 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
90 : real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
91 : real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
92 : real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
93 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
94 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
95 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
96 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
97 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
98 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
99 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
100 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
101 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
102 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
103 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
104 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
105 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
106 : parameter(alpha = 0.1875D0, beta = 0.125D0)
107 : c-----------------------
108 : gm1=ga-1.0d0
109 : c-----------------------
110 : n_x=nvect(1)
111 : n_y=nvect(2)
112 : t_x=tvect(1)
113 : t_y=tvect(2)
114 : c----------------------------
115 : c----------------------------
116 : rold_l=wvec_l(1)
117 : uold_l=wvec_l(2)
118 : vold_l=wvec_l(3)
119 : pold_l=wvec_l(4)
120 : c-----------------------
121 : rold_r=wvec_r(1)
122 : uold_r=wvec_r(2)
123 : vold_r=wvec_r(3)
124 : pold_r=wvec_r(4)
125 : c------------------------------------------------------------------
126 : c Computation of the specific total energy on the left and right.
127 : c------------------------------------------------------------------
128 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
129 : eold_l=eold_l+pold_l/(gm1*rold_l)
130 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
131 : eold_r=eold_r+pold_r/(gm1*rold_r)
132 : c-------------------------------------------------------------------
133 : c Computation of the speed of sound and its derivatives;
134 : c numerical speed of sound at the interface is taken as an average
135 : c of the speeds of sounds of the neighbouring cells
136 : c-------------------------------------------------------------------
137 : aleft=dsqrt(ga*pold_l/rold_l)
138 : arigh=dsqrt(ga*pold_r/rold_r)
139 : am=0.5d0*(aleft+arigh)
140 : c-------------------------------------------------------------------
141 : damr_r=-arigh/(4.0d0*rold_r)
142 : damu_r=0.0d0
143 : damv_r=0.0d0
144 : damp_r=ga/(4.0d0*arigh*rold_r)
145 : c-----------------------
146 : damr_l=-aleft/(4.0d0*rold_l)
147 : damu_l=0.0d0
148 : damv_l=0.0d0
149 : damp_l=ga/(4.0d0*aleft*rold_l)
150 : c-------------------------------------------------------------------
151 : c Computing numerical Mach number and its derivatives,
152 : c see p.370, under (A1).
153 : c-------------------------------------------------------------------
154 : un_l=uold_l*n_x+vold_l*n_y
155 : un_r=uold_r*n_x+vold_r*n_y
156 : ut_l=uold_l*t_x+vold_l*t_y
157 : ut_r=uold_r*t_x+vold_r*t_y
158 : c-------------------------------------------------------------------
159 : ml=un_l/am
160 : mr=un_r/am
161 : c-------------------------------------------------------------------
161 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
162 : c see (19a) and (19b), p.367.
163 : c-------------------------------------------------------------------
164 : if(dabs(ml) .ge. 1.0d0) then
165 : Mplus=(ml+dabs(ml))/2.0d0
166 : else
167 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
168 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
169 : endif
170 : c-------------------------------------------------------------------
171 : if(dabs(mr) .ge. 1.0d0) then
172 : Mmin=(mr-dabs(mr))/2.0d0
173 : else
174 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
175 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
176 : endif
177 : c-------------------------------------------------------------------
178 : c Derivatives of ml and mr with respect to both sets of primitive
179 : c variables: left and right.
180 : c-------------------------------------------------------------------
181 : temp_l=-un_l/(am*am)
182 : temp_r=-un_r/(am*am)
183 : c--------
184 : dmlr_l=temp_l*damr_l
185 : dmlr_r=temp_l*damr_r
186 : dmrr_l=temp_r*damr_l
187 : dmrr_r=temp_r*damr_r
188 : c--------
189 : dmlu_l=n_x/am+temp_l*damu_l
190 : dmlu_r=temp_l*damu_r
191 : dmru_l=temp_r*damu_l
192 : dmru_r=n_x/am+temp_r*damu_r
193 : c--------
194 : dmlv_l=n_y/am+temp_l*damv_l
195 : dmlv_r=temp_l*damv_r
196 : dmrv_l=temp_r*damv_l
197 : dmrv_r=n_y/am+temp_r*damv_r
198 : c--------
199 : dmlp_l=temp_l*damp_l
200 : dmlp_r=temp_l*damp_r
201 : dmrp_l=temp_r*damp_l
202 : dmrp_r=temp_r*damp_r
203 : c-----------------------------------------------------------
204 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
205 : c-----------------------------------------------------------
206 : mmid=Mplus+Mmin
207 : c-----------------------------------------------------------
208 : c Computing the derivatives of M+ and M-
209 : c-----------------------------------------------------------
210 : if(ml .ge. 1.0d0) then
211 : dMpr_l=dmlr_l
212 : dMpu_l=dmlu_l
213 : dMpv_l=dmlv_l
214 : dMpp_l=dmlp_l
215 : c--------------------
216 : dMpr_r=dmlr_r
217 : dMpu_r=dmlu_r
218 : dMpv_r=dmlv_r
219 : dMpp_r=dmlp_r
220 : else
221 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
222 : temph=(ml+1.0d0)/2.0d0
223 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
224 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
225 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
226 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
227 : c--------------------
228 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
229 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
230 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
231 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
232 : else
233 : dMpr_l=0.0d0
234 : dMpu_l=0.0d0
235 : dMpv_l=0.0d0
236 : dMpp_l=0.0d0
237 : c---------------------
238 : dMpr_r=0.0d0
239 : dMpu_r=0.0d0
240 : dMpv_r=0.0d0
241 : dMpp_r=0.0d0
242 : endif
243 : endif
244 : c-----------------------------------------------------------
245 : if(mr .ge. 1.0d0) then
246 : dMmr_l=0.0d0
247 : dMmu_l=0.0d0
248 : dMmv_l=0.0d0
249 : dMmp_l=0.0d0
250 : c---------------------
251 : dMmr_r=0.0d0
252 : dMmu_r=0.0d0
253 : dMmv_r=0.0d0
254 : dMmp_r=0.0d0
255 : else
256 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
257 : temph=(-mr+1.0d0)/2.0d0
258 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
259 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
260 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
261 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
262 : c--------------------
263 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
264 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
265 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
266 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
267 : else
268 : dMmr_l=dmrr_l
269 : dMmu_l=dmru_l
270 : dMmv_l=dmrv_l
271 : dMmp_l=dmrp_l
272 : c--------------------
273 : dMmr_r=dmrr_r
274 : dMmu_r=dmru_r
275 : dMmv_r=dmrv_r
276 : dMmp_r=dmrp_r
277 : endif
278 : endif
279 : c-----------------------------------------------------------------
280 : c computing the derivatives of m_{1/2} (notation as in the paper)
281 : c-----------------------------------------------------------------
282 : dmir_l=dMpr_l+dMmr_l
283 : dmir_r=dMpr_r+dMmr_r
284 : c-------------
285 : dmiu_l=dMpu_l+dMmu_l
286 : dmiu_r=dMpu_r+dMmu_r
287 : c-------------
288 : dmiv_l=dMpv_l+dMmv_l
289 : dmiv_r=dMpv_r+dMmv_r
290 : c-------------
291 : dmip_l=dMpp_l+dMmp_l
292 : dmip_r=dMpp_r+dMmp_r
293 : c----------------------------------------------------------------
294 : c computing the main convective variables and their derivatives
295 : c mpl_m is m^{+}_{1/2} (paper's notation) and
296 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
297 : c----------------------------------------------------------------
298 : if(mmid .ge. 0.0d0) then
299 : mpl_m = mmid
300 : d2mr_l=dmir_l
301 : d2mu_l=dmiu_l
302 : d2mv_l=dmiv_l
303 : d2mp_l=dmip_l
304 : c------------
305 : d2mr_r=dmir_r
306 : d2mu_r=dmiu_r
307 : d2mv_r=dmiv_r
308 : d2mp_r=dmip_r
309 : c------------
310 : else
311 : mpl_m = 0.0d0
312 : d2mr_l=0.0d0
313 : d2mu_l=0.0d0
314 : d2mv_l=0.0d0
315 : d2mp_l=0.0d0
316 : c------------
317 : d2mr_r=0.0d0
318 : d2mu_r=0.0d0
319 : d2mv_r=0.0d0
320 : d2mp_r=0.0d0
321 : endif
322 : c---------------------------------------------
323 : if(mmid .le. 0.0d0) then
324 : mmin_m = mmid
325 : d3mr_l=dmir_l
326 : d3mu_l=dmiu_l
327 : d3mv_l=dmiv_l
328 : d3mp_l=dmip_l
329 : c------------
330 : d3mr_r=dmir_r
331 : d3mu_r=dmiu_r
332 : d3mv_r=dmiv_r
333 : d3mp_r=dmip_r
334 : c------------
335 : else
336 : mmin_m = 0.0d0
337 : d3mr_l=0.0d0
338 : d3mu_l=0.0d0
339 : d3mv_l=0.0d0
340 : d3mp_l=0.0d0
341 : c------------
342 : d3mr_r=0.0d0
343 : d3mu_r=0.0d0
344 : d3mv_r=0.0d0
345 : d3mp_r=0.0d0
346 : endif
347 : c---------------------------------------------------------------
347 : c Computing the calligraphic P+ and P- with their derivatives,
348 : c see (21a) & (21b) on p.368.
349 : c---------------------------------------------------------------
350 : if(ml .ge. 1.0d0) then
351 : Pplus = 1.0d0
352 : else
353 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
354 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
355 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
356 : else
357 : Pplus = 0.0d0
358 : endif
359 : endif
360 : c---------------------------------------------------------------
361 : if(mr .ge. 1.0d0) then
362 : Pmin = 0.0d0
363 : else
364 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
365 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
366 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
367 : else
368 : Pmin = 1.0d0
369 : endif
370 : endif
371 : c---------------------------------------------------------------
372 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
373 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
374 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
375 : c--------------
376 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
377 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
378 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
379 : c---------------------------------------------------------------
380 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
381 : dPpr_l=brac_l*dmlr_l
382 : dPpr_r=brac_l*dmlr_r
383 : c------------
384 : dPpu_l=brac_l*dmlu_l
385 : dPpu_r=brac_l*dmlu_r
386 : c------------
387 : dPpv_l=brac_l*dmlv_l
388 : dPpv_r=brac_l*dmlv_r
389 : c------------
390 : dPpp_l=brac_l*dmlp_l
391 : dPpp_r=brac_l*dmlp_r
392 : c------------
393 : else
394 : dPpr_l=0.0d0
395 : dPpr_r=0.0d0
396 : c-----------
397 : dPpu_l=0.0d0
398 : dPpu_r=0.0d0
399 : c-----------
400 : dPpv_l=0.0d0
401 : dPpv_r=0.0d0
402 : c-----------
403 : dPpp_l=0.0d0
404 : dPpp_r=0.0d0
405 : c-----------
406 : endif
407 : c---------------------------------------------------------------
408 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
409 : dPmr_l=brac_r*dmrr_l
410 : dPmr_r=brac_r*dmrr_r
411 : c------------
412 : dPmu_l=brac_r*dmru_l
413 : dPmu_r=brac_r*dmru_r
414 : c------------
415 : dPmv_l=brac_r*dmrv_l
416 : dPmv_r=brac_r*dmrv_r
417 : c------------
418 : dPmp_l=brac_r*dmrp_l
419 : dPmp_r=brac_r*dmrp_r
420 : c------------
421 : else
422 : dPmr_l=0.0d0
423 : dPmr_r=0.0d0
424 : c-----------
425 : dPmu_l=0.0d0
426 : dPmu_r=0.0d0
427 : c-----------
428 : dPmv_l=0.0d0
429 : dPmv_r=0.0d0
430 : c-----------
431 : dPmp_l=0.0d0
432 : dPmp_r=0.0d0
433 : c-----------
434 : endif
435 : c-------------------------------------------------------------------
436 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
437 : c-------------------------------------------------------------------
438 : pmid=Pplus*pold_l + Pmin*pold_r
439 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
440 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
441 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
442 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
443 : c----------------------------
444 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
445 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
446 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
447 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
448 : c---------------------------------------------------------------------
449 : c computing JACOBIAN as a derivative of the numerical flux function
450 : c with respect to the primitive variables.
451 : c Notation: jl(i,j) --- is the derivative of the i-component of the
452 : c flux function with respect to the j-component of the
453 : c vector of primitive variables of the left state.
454 : c jr(i,j) --- is the derivative of the i-component of the
455 : c flux function with respect to the j-component of the
456 : c vector of primitive variables of the right state.
457 : c---------------------------------------------------------------------
458 : f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
459 : c---------------------------------------------------------------------
460 : jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
461 : jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
462 : jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
463 : jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
464 : jl(1,4)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
465 : c------------------------------------
466 : jr(1,1)=damr_r*f(1)/am+am*(d2mr_r*rold_l+mmin_m)
467 : jr(1,1)=jr(1,1)+am*d3mr_r*rold_r
468 : jr(1,2)=damu_r*f(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
469 : jr(1,3)=damv_r*f(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
470 : jr(1,4)=damp_r*f(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
471 : c------------------------------------
472 : c------------------------------------
473 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
474 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
475 : f(2)=n_x*(am*br1+pmid)+am*t_x*br2
476 : c------------------
477 : det=n_x*t_y-n_y*t_x
478 : c---------------------------------------------------------
479 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
480 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
481 : jl(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
482 : jl(2,1)=jl(2,1)+t_x*damr_l*br2+am*t_x*br4
483 : c-------------------
484 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
485 : br3=br3+d3mu_l*rold_r*un_r
486 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y)/det
487 : br4=br4+d3mu_l*rold_r*ut_r
488 : jl(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
489 : jl(2,2)=jl(2,2)+t_x*damu_l*br2+am*t_x*br4
490 : c-------------------
491 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x)/det
492 : br3=br3+d3mv_l*rold_r*un_r
493 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
494 : br4=br4+d3mv_l*rold_r*ut_r
495 : jl(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
496 : jl(2,3)=jl(2,3)+t_x*damv_l*br2+am*t_x*br4
497 : c-------------------
498 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
499 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
500 : jl(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
501 : jl(2,4)=jl(2,4)+t_x*damp_l*br2+am*t_x*br4
502 : c-------------------------------------------------------------
503 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
504 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
505 : jr(2,1)=n_x*(damr_r*br1+am*br3+dpir_r)
506 : jr(2,1)=jr(2,1)+t_x*damr_r*br2+am*t_x*br4
507 : c-------------------
508 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
509 : br3=br3+d3mu_r*rold_r*un_r
510 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
511 : br4=br4+d3mu_r*rold_r*ut_r
512 : jr(2,2)=n_x*(damu_r*br1+am*br3+dpiu_r)
513 : jr(2,2)=jr(2,2)+t_x*damu_r*br2+am*t_x*br4
514 : c-------------------
515 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
516 : br3=br3+d3mv_r*rold_r*un_r
517 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
518 : br4=br4+d3mv_r*rold_r*ut_r
519 : jr(2,3)=n_x*(damv_r*br1+am*br3+dpiv_r)
520 : jr(2,3)=jr(2,3)+t_x*damv_r*br2+am*t_x*br4
521 : c-------------------
522 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
523 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
524 : jr(2,4)=n_x*(damp_r*br1+am*br3+dpip_r)
525 : jr(2,4)=jr(2,4)+t_x*damp_r*br2+am*t_x*br4
526 : c-------------------------------------------------------------
527 : c-------------------------------------------------------------
528 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
529 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
530 : f(3)=n_y*(am*br1+pmid)+am*t_y*br2
531 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
532 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
533 : jl(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
534 : jl(3,1)=jl(3,1)+t_y*damr_l*br2+am*t_y*br4
535 : c-------------------
536 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
537 : br3=br3+d3mu_l*rold_r*un_r
538 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y/det)
539 : br4=br4+d3mu_l*rold_r*ut_r
540 : jl(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
541 : jl(3,2)=jl(3,2)+t_y*damu_l*br2+am*t_y*br4
542 : c-------------------
543 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x/det)
544 : br3=br3+d3mv_l*rold_r*un_r
545 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
546 : br4=br4+d3mv_l*rold_r*ut_r
547 : jl(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
548 : jl(3,3)=jl(3,3)+t_y*damv_l*br2+am*t_y*br4
549 : c-------------------
550 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
551 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
552 : jl(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
553 : jl(3,4)=jl(3,4)+t_y*damp_l*br2+am*t_y*br4
554 : c-------------------------------------------------------------
555 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
556 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
557 : jr(3,1)=n_y*(damr_r*br1+am*br3+dpir_r)
558 : jr(3,1)=jr(3,1)+t_y*damr_r*br2+am*t_y*br4
559 : c-------------------
560 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
561 : br3=br3+d3mu_r*rold_r*un_r
562 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
563 : br4=br4+d3mu_r*rold_r*ut_r
564 : jr(3,2)=n_y*(damu_r*br1+am*br3+dpiu_r)
565 : jr(3,2)=jr(3,2)+t_y*damu_r*br2+am*t_y*br4
566 : c-------------------
567 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
568 : br3=br3+d3mv_r*rold_r*un_r
569 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
570 : br4=br4+d3mv_r*rold_r*ut_r
571 : jr(3,3)=n_y*(damv_r*br1+am*br3+dpiv_r)
572 : jr(3,3)=jr(3,3)+t_y*damv_r*br2+am*t_y*br4
573 : c-------------------
574 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
575 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
576 : jr(3,4)=n_y*(damp_r*br1+am*br3+dpip_r)
577 : jr(3,4)=jr(3,4)+t_y*damp_r*br2+am*t_y*br4
578 : c-------------------------------------------------------------
579 : c ------ f44444444444444444444444444444 ---------
580 : c-------------------------------------------------------------
581 : hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1
582 : hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1
583 : f(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
584 : c---------------------
585 : br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
586 : br1=br1+d3mr_l*hr_r
587 : jl(4,1)=damr_l*f(4)/am+am*br1
588 : c---------------------
589 : br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
590 : br1=br1+d3mu_l*hr_r
591 : jl(4,2)=damu_l*f(4)/am+am*br1
592 : c---------------------
593 : br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
594 : br1=br1+d3mv_l*hr_r
595 : jl(4,3)=damv_l*f(4)/am+am*br1
596 : c---------------------
597 : br1=d2mp_l*hr_l+mpl_m*ga/gm1
598 : br1=br1+d3mp_l*hr_r
599 : jl(4,4)=damp_l*f(4)/am+am*br1
600 : c----------------------------------------------------------
601 : c----------------------------------------------------------
602 : br1=d2mr_r*hr_l+mmin_m*(uold_r*uold_r+vold_r*vold_r)/2.0d0
603 : br1=br1+d3mr_r*hr_r
604 : jr(4,1)=damr_r*f(4)/am+am*br1
605 : c---------------------
606 : br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
607 : br1=br1+d3mu_r*hr_r
608 : jr(4,2)=damu_r*f(4)/am+am*br1
609 : c---------------------
610 : br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
611 : br1=br1+d3mv_r*hr_r
612 : jr(4,3)=damv_r*f(4)/am+am*br1
613 : c---------------------
614 : br1=d2mp_r*hr_l+mmin_m*ga/gm1
615 : br1=br1+d3mp_r*hr_r
616 : jr(4,4)=damp_r*f(4)/am+am*br1
617 : c-------------------------------------------------------------
618 : c-------------------------------------------------------------
619 : return
620 : end
621 :
622 :
623 :
624 :
625 :
626 :
627 :
628 :
629 :
630 :
631 :
632 :
633 :
634 :
635 :
636 :
637 :
638 :
639 :
640 :
641 :
© Cast3M 2003 - All rights reserved.
Disclaimer