so here is my code. Clearly it doesn't work right now, and it's very long. Getting it to calculate velocities would be a good start, and i'm pretty sure my flagging is not correct.
$xlength=10; $ylength=4; $N=100; my @u; #Horizontal Velocity my @v; #Vertical Velocity my @p; my @F; my @G; my $eps=.001; my @grid; $imax=25; $jmax=10; for $i (0..$imax) { for $j (0..$jmax) { $u[$i][$j]=0; $v[$i][$j]=0; $p[$i][$j]=0; }} for $i (1..$imax){ for $j (1..$jmax) { $u[0][$j]=0; $v[$i][$jmax]=0; $v[0][$j]=$v[1][$j]; $u[$i][$jmax+1]=$u[$i][$jmax]; }} for $i (1..$imax) { for $j (1..$jmax) { $u[$i][1]=100; $u[$imax][$j]=100; $v[$i][1]=100; $v[$imax][$j]=100; $p[$i][1]=100; $p[$imax][$j]=100; }} $gx=0; $gy=-1; $delx=1; $dely=1; $RE=10; $time_step=.5; print "$v[10][1]\n"; ##### Assigning Flags###### #Begin Time loop#### for ($t=0; $t <$tend; $t +=$time_step) { for $i (0..$N) { for $j (0..$N) { if ($v[$i][$j] and $u[$i][$j] > 0 ){$flag=CF;} if ($v[$i][$j] and $u[$i][$j] = 0) {$flag=CE;} for $flag (CF) { if( $v[$i][$j-1]=0 or $u[$i][$j-1]=0){$flag='C_N';} if( $v[$i][$j+1]=0 or $u[$i][$j+1]=0) {$flag='C_S';} if ($u[$i-1][$j]=0 or $v[$i-1][$j]=0){$flag='C_W';} if ($u[$i+1][$j]=0 or $v[$i+1][$j]=0) {$flag='C_E';} if ($u[$i][$j-1]=0 and $u[$i][$j+1]=0 or $v[$i][$j-1]=0 and $v[$i][$ +j+1]=0){$flag='C_NS';} if ($u[$i-1][$j]=0and $u[$i+1][$j]=0 or $v[$i-1][$j]=0and $v[$i+1][$j] +=0){$flag='C_WE';} if ($u[$i][$j-1]=0 and $u[$i-1][$j]=0 or $v[$i][$j-1]=0 and $v[$i-1][$ +j]=0){$flag='C_NW';} if ($u[$i][$j-1]=0 and $u[$i+1][$j]=0 or $v[$i][$j-1]=0 and $v[$i+1][ +$j]=0){$flag='C_NE';} if ($u[$i][$j+1]=0 and $u[$i-1][$j]=0 or $v[$i][$j+1]=0 and $v[$i-1][ +$j]=0){$flag='C_SW';} if ($u[$i][$j+1]=0 and $u[$i+1][$j]=0 or $v[$i][$j+1]=0 and $v[$i+1][ +$j]=0){$flag='C_SE';} if ($u[$i][$j-1]=0 and $u[$i-1][$j]=0 and $u[$i+1][$j]=0 or $v[$i][$j +-1]=0 and $v[$i-1][$j]=0 and $v[$i+1][$j]=0){$flag='C_NWE';} if ($u[$i][$j-1]=0 and $u[$i][$j+1]=0 and $u[$i+1][$j]=0 or $v[$i][$j +-1]=0 and $v[$i][$j+1]=0 and $v[$i+1][$j]=0){$flag='C_NSE';} if ($u[$i][$j-1]=0 and $u[$i][$j+1]=0 and $u[$i-1][$j]=0 or $v[$i][$j +-1]=0 and $v[$i][$j+1]=0 and $v[$i-1][$j]=0){$flag='C_NSW';} if ($u[$i+1][$j]=0 and $u[$i][$j+1]=0 and $u[$i-1][$j]=0 or $v[$i+1][$ +j]=0 and $v[$i][$j+1]=0 and $v[$i-1][$j]=0){$flag='C_SWE';} if ($u[$i][$j-1]=0 and $u[$i-1][$j]=0 and $u[$i][$j+1]=0 and $u[$i+1] +[$j]=0 or $v[$i][$j-1]=0 and $v[$i-1][$j]=0 and $v[$i][$j+1]=0 and $ +v[$i+1][$j]=0){$flag='C_NWSE';} for $i (1..$imax) { for $j (1..$jmax) { if ($flag eq 'C_N', 'C_S', 'C_E', 'C_W'){ $p[$i][$j]=2/$RE*(($u[$i][$j]-$u[$i-1][$j])/$delx); $u[$i][$j]=$u[$i-1][$j]-($delx/$dely)*($v[$i][$j]-$v[$i][$j-1]); $v[$i+1][$j-1]=$v[$i][$j-1]-($delx/$dely)*($u[$i][$j]-$u[$i][$j-1]); print "$u[10][1]\n" } if ($flag='C_NE','C_SW'){ $u[$i][$j]=$u[$i-1][$j]; $v[$i][$j-1]=$v[$i][$j]; $p[$i][$j]= 1/2*$RE*((($u[$i][$j+1]+$u[$i-1][$j+1]-$u[$i][$j]-$u[$i-1] +[$j])/$dely)+(($v[$i][$j]+$v[$i][$j-1]-$v[$i-1][$j]-$v[$i-1][$j-1])/$ +delx)); } if ($flag='C_NW','C_SE') { $u[$i][$j]=$u[$i+1][$j]; $v[$i][$j+1]=$v[$i][$j]; $p[$i][$j]= -1/2*$RE*(($u[$i][$j+1]+$u[$i-1][$j+1]-$u[$i][$j]-$u[$i-1] +[$j])/$dely)+(($v[$i][$j]+$v[$i][$j-1]-$v[$i-1][$j]-$v[$i-1][$j-1])/$ +delx); } if ($flag='C_WE', 'C_NS') { $u[$i][$j]=$u[$i][$j]+$time_step*$gx; $u[$i-1][$j]=$u[$i-1][$j]+$time_step*$gx; $v[$i][$j]=$v[$i][$j]+$time_step*$gy; $v[$i][$j-1]=$v[$i][$j-1]+$time_step*$gy; $p[$i][$j]=0; } if ($flag='C_NSE','C_NSW','C_SWE','C_NWE'){ $p[$i][$j]=0; } if ($flag='C_NWSE'){ $u[$i][$j]=$u[$i][$j+1]; $p[$i][$j]=0; $v[$i][$j]=$v[$i-1][$j]; } } } for $i (1..$imax-1){ for $j(1..$jmax){ $F[$n][$i][$j]=$u[$i][$j]+$time_step*(1/$RE*(($u[$i][$j]/($xlength*$xl +ength))+($u[$i][$j]/($ylength*$ylength)))-(($u[$i][$j]*$u[$i][$j])/$x +length)-(($u[$i][$j]*$v[$i][$j])/$ylength)+$gx); $G[$n][$i][$j]=$v[$i][$j]+$time_step*(1/$RE*(($u[$i][$j]/($xlength*$xl +ength))+($u[$i][$j]/($ylength*$ylength)))-(($u[$i][$j]*$v[$i][$j])/$x +length)-(($u[$i][$j]*$u[$i][$j])/$ylength)+$gy); } } ##### more boundary conditions for RHS ###### for $i (1..$imax) { for $j (1..$jmax){ $p[0][$j]=$p[1][$j]; $p[$imax+1][$j]=$p[$imax][$j]; $p[$i][0]=$p[$i][1]; $p[$i][$jmax+1]=$p[$i][$jmax]; $Fn[0][$j]=$u[0][$j]; $Fn[$imax][$j]=$u[$imax][$j]; $Gn[$i][0]=$v[$i][0]; $Gn[$i][$jmax]=$v[$i][$jmax]; }} ####RHS Calc### for $i (1..$imax){ for $j(1..$jmax){ $RHS[$i][$j]=(1/$time_step)*(($Fn[$i][$j]-$Fn[$i-1][$j])/$delx+($Gn[$i +][$j]-$Gn[$i][$j-1])/$dely); if ($i=1){ $eW=0; }elsif ($i>1) { $eW=1; }elsif ($i<$imax){ $eE=1; }elsif ($i=$imax) { $eE=0; } if ($j=1){ $eS=0; }elsif ($j>1){ $eS=1; }elsif ($j<$jmax){ $eN=1; }elsif ($j=$jmax){ $eN=0; } } } ####SOR ##### #$itmax=20; #$omg=1; #do #{ #for $it (1..$itmax){ #for $i (1..$imax) { #for $j (1..$jmax) { #$P[$i][$j]=(1-$omg)*$p[$i][$j]+($omg/(($eE+$eW)/($delx*$delx))+($eN+$ +eS)/($dely*$dely)))*((($eE*$p[$i+1][$j]+$eW*$P[$i-1][$j])/($delx*$del +x))+(($eN*$p[$i][$j+1]+$eS*$P[$i][$j-1])/($dely*$dely)-$RHS[$i][$j]); #$rit[$i][$j]=((($eE*($p[$i+1][$j]-$p[$i][$j])-$eW*($p[$i][$j]-$p[$i-1 +][$j]))/($delx*$delx))+(($eN*($p[$i][$j+1]-$p[$i][$j])-$eS*($p[$i][$j +]-$p[$i][$j-1])/($dely*$dely))-$RHS[$i][$j])); #}while ($it<$itmax or abs($rit) gt $eps); #} #} #} }}} }
In reply to Re^4: Particle movement question
by koda123
in thread Particle movement question
by koda123
| For: | Use: | ||
| & | & | ||
| < | < | ||
| > | > | ||
| [ | [ | ||
| ] | ] |